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ABSTRACT 



This paper investigates some of the difficulties 
commonly encountered in the application of pattern recog- 
nition techniques to the identification of non-cooperative 
signals. Practical problems of feature selection and analy- 
sis are explored for a data base of actual bauded signal 
transition times. Analysis of a set of features based on 
the raster display of these transition times suggests that 
the use of fast Fourier transform techniques may yield 
useful identifying parameters not exploitable by the current 
process of visual comparison of signal raster patterns. 
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I. INTRODUCTION 



The manual measurement and analysis of signal parameters 
leading to identification of the signal source can require 
more time and manpower than results justify. Such a process 
may be made less costly by the application of computer tech- 
nology to measure the traditional parameters in a manner 
relatively free of individual bias and to aid in the analy- 
sis and identification procedures. An additional result of 
the computer application is the capability to calculate 
parameters unmeasurable by manual techniques. This thesis 
investigates the process by which such parameters might be 
obtained and used in an existing pattern recognition scheme. 
The discussion progresses from the general principles of 
Bayesian classification to the problem of finding features 
which are useful in classifying signals from a specific 
data base. 

Section II of this thesis deals with the theory of 
pattern recognition which underlies most automated classi- 
fication schemes. Section III discusses some of the limi- 
tation imposed oxi this theory in a practical signal identi- 
fication problem. Section IV is devoted to the problems 
v/hich arise in the attempt to choose a set of features useful 
and meaningful to the identification process. Section V 
addressed the particular problems of choosing a set of 
useful features based on the raster scan display of bauded 
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signal transition times. Included in this section are 
the motivation for the initial choice of parameters and 
the techniques used for both measurement and analysis. The 
conclusions reached by this investigator, as well as those 
areas deserving further study, are summarized in Section VI. 
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II. BAYESIAN PATTERN CLASSIFICATION 



The process of mathematical pattern recognition forma- 
lizes the methods by which the measured parameters or fea- 
tures of a group of samples may be used to identify those 
samples as belonging to one of several general classes. 

Such a process may be considered to be composed of two 
parts. The first part is a learning procedure in v;hich the 
general "pattern" of each of the various classes is des- 
cribed in terms of the similarities and differences in the 
sample features. The second part is a decision process in 
v;hich a feature-based algorithm is developed from the pattern 
description and applied to the identification or classifi- 
cation of nev; samples. Both parts of the process are based 
on the theory of statistical inference and lead to a pro- 
babilistic interpretation of the class membership. Closely 
related to such an interpretation is the Bayesian model of 
decision theory. 

In the Bayesian model, the first step of a pattern recog- 
nition scheme, the learning process, is that technique by 
which a-priori class-conditional feature probabilities are 
determined. The K^-individual d-dimensional feature vectors 
^2^^' ••• / distinct classes are 

used to estimate m class conditional probability density 
functions on the feature space 
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• • • 



(II-l) 



p(x|c^) ~ P(2il2ii 



(i) 



1 



Knowing these distributions for each class and the overall 
class probabilities, p(C^), one may obtain the feature 
conditional class probabilities: 



p(C^ lx) 



P(C. )p(xlc.) 
p(x) “ 



p(C^)p(x|C^) 

__ 

Z p(c . )p(x|c. ) 

i=l ^ ^ 



(II-2) 



Assuming an equal cost of misclassification , the Bayesian 
decision rule chooses that class V7hich maximizes the feature 
conditional probabilities of Equation II"2. That is, if 



p(C^|x) >_ p(Cj|x) for all j (II-3) 

the sample characterized by the feature vector x is assigned 

to class C . . 

± 

The above technique generates a series of decision 
regions in the d-dimensional feature space. Once the a-priori 
probabilities are determined, the class selection is simply 
an m-valued function of the feature vector. This is essen- 
tially the goal of the pattern recognition process. 

The probabilistic Bayesian decision process discussed 
above is mathematically elegant and optimal in the sense 
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that if the a-priori probabilities are known, it minimizes 
the expected risk of misclassification . Unfortunately, 
its implementation suffers major difficulties in both the 
learning and the decision phases of a statistical pattern 
recognition process. 

A. PATTERN LEARNING 

In the learning portion of a pattern recognition problem, 
one normally deals with a finite number of samples which 
constitute a training set. These samples are the basis from 
which the investigator is obliged to construct the necessary 
a-priori class-conditional probability distribution functions. 
Thus the well-defined probability structure assumed by 
Bayes* decision rule is in reality a statistical inference 
based on the data obtained from the training set and the 
investigator’s best guess of the functional form of this 
data. Dependent on hov; m.uch of the general form of the 
probability density function is presupposed, the estimation 
techniques are characterized as parairvetric or non-parametric . 

Typically, parametric estimation assumes that the features 
of the training samples follow one of the well-known proba- 
bility distributions. The samples of the training set are 
used to estimate only a few essential parameters of this 
distribution. For example, one might use the training sam- 
ples to estimate the means and variances for an assumed 
multivariate normal distribution. The non-parametric tech- 
nique on the other hand seeks to estimate directly the 
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probability density function as a superposition of "potential 
function" terras arising from each of the training set samples. 

The accuracy of the inferred density function is strongly 
dependent on the type of estimation technique used, and may 
be severely limited by the size of the training set. The 
applicability and limitations of both non-parametric and 
parametric estimation are discussed further in Appendix A. 
Regardless of the type of estimation, the process involved 
is an example of supervised learning. It presupposes a 
set of correctly labeled samples from which the estimates 
may be constructed. Information external to the features 
being measured has been used to identify the members of 
each class. 

B. DECISION RULES 

The previous discussion of the Bayesian approach to 
signal classification already suggests one method of ob- 
taining a decision rule. For each new (unclassified) 
feature vector use the a-priori probability inferred from 
the taining set to calculate the feature-conditional pro- 
bability of its membership in each class; then choose that 
class for which this probability is greatest. 

Such a maximal likelihood criteria is intuitively 
satisfying and can be formulated to lead to classifications 
which are optimal in the sense of minimizing the average 
risk of misclassification . Unfortunately a naive implemen- 
tation of such a technique leads to computational inefficiencies 
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which often render the process impractical in terms of 
memory and time constraints. 

More frequently the decision rule chosen in a practical 
pattern recognition problem is suboptimal in the Bayesian 
sense, providing a method of partitioning the feature-space 
of the samples without resorting to the calculation of the 
maximal likelihood criteria. These suboptimal schemes 
use simplifying assumptions regarding the nature of the 
underlying class-conditional probability distributions to 
arrive at class boundaries which are relatively simple func- 
tions of the features themselves. If there is little overlap 
in the various class-conditional probability density func- 
tions, the increase in the expected misclassification rate 
as a result of the stiboptimal decision is minimal. 

For example, consider the case of the one-dimensional 
distributions for the two equiprobable classes shown in 
Figures la and lb. Suppose that some suboptimal decision 
rule results in decision boundary 

B’ = B . - e (II- 4) 

opt 

The shaded areas of Figures la and lb indicate the resultant 
misclassification rates. 

If there is substantial separation of the class distri- 
butions, the differential error in classification as a result 
of boundary misplacement, area D, is negligible. In general. 



17 




I 





FIGURE 1. Misclassification Rates Due ot Suboptinal 
Decision Boundaries 



18 



with decreasing overlap of unimodal class distributions, 
not only does the total misclassification error decrease, 
but so does its sensitivity to the placement of the decision 
boundary. As a practical matter, the distinction between 
a five percent and a two percent misclassification rate 
seems unneccessarily fine if the probability distributions 
were originally inferred from ten samples each. 

In many cases the sub-optimal classification schemes 
follow rather naturally from the method used to infer the 
a-priori distributions. Unimodal feature distributions 
give rise to classification schemes based on a test sample's 
minimal "distance" from estimated class means. More complex 
feature distributions involve the use of a distance measure 
in the feature space to establish the K training set members 
closest to the test sample; but classification based on the 
weighted vote of these K-nearest neighbors, seems closer in 
spirit to non-parametric probability estimation. Both 
these classification methods are explored more fully in 
Appendix B . 
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III. PRACTICAL CONSIDERATIONS IN THE IDENTIFICATION 
OF NON-COOPERATIVE SIGNALS 



The foregoing discussion of the basic techniques of 
mathematical pattern recognition is intended to serve as 
an outline of approaches to the general problems of classi- 
fication. Identification of the originator of a bauded 
signal, while bearing considerable similarity to Bayesian 
classification, is complicated by a number. of factors not 
normally arising in such problems as the identification of 
handvr?ritten characters . 

In many of the classical problems of mathematical pat- 
tern recognition, samples are obtained from a statistical 
universe v;hich is well-defined in the sense that the number 
of classes to which the samples may belong is finite and 
previoiisly ascertained. Additionally, the parameters or 
features used to typify these classes are normally time- 
invariant. The set used for training is accurately classi- 
fied and sufficiently large that meaningful statistical 
estimates of the parameters and their distributions may 
be obtained. 

Under these conditions, the processes of Section I 
find their greatest application. Unfortunately, they are 
seldom encountered in identifying the signals of a non- 
cooperative originator. The data base used for a training 
set does not normally provide arbitrarily accurate ground 
truth classification. In practice, even the number of 
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actual originator classes may not be known. In many cases, 
the paucity of signal intercept precludes more than the 
grossest statistical parameter estimation. The parameters 
themselves are time-dependent since they are affected by 
component aging, time-dependent signal propagation condi- 
tions, and unknown external maintenance and adjustment. 

The advantages of long-term averaging to obtain statistical 
convergence, i.e. accurate a-priori estimates, may be 
obscured by short-term perturbations for which adequate 
correction may be prohibitively costly or undefinable. 

In summary, the investigator proposing to use pattern 
recognition techniques to identify the originator of hostile 
signals is faced with a series of constraints involving 
small sample size, time-varying signcil parameters and class 
probabilities, inaccurate ground truth, and an indeterminate 
number of classes with which signals may be associated. 

Such constraints impose severe limitations upon the pattern 
recognition techniques V7hich may be successfully applied to 
the problem. The following sections indicate some of the 
procedures which deal with these limitations . 

A. SMALL SAMPLE SIZE 

The difficulties arising from small sample size are most 
immediately apparent when one attempts to increase the number 
of features in order to facilitate the separation of classes. 
The dominant effect is one of inaccuracy arising from under- 
sampling the required d-dimensional a-priori probability 
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distributions. V'Jhen this undersampling occurs, the most 
likely result is that increasing the number of features be- 
comes counterproductive to "good" classification. In the 
practical multi-category problem, it is not unusual to 
encounter an initial set as many as fifty features which 
the designer believes to be useful in identification. The 
number of taining samples is often too small to allow 
meaningfiil Parzen window or other non-pararaetric density 
estimation techniques over this number of dimensions. 

An alternative procedure of characterizing the a-priori 
distribution function by use of the most general multi- 
variate normal form requires estimation of a "d x d" non- 
singular covariance matrix. This imposes an algebraic re- 
quirement of at least d+1 samples. However, experience in 
the practical estimation of covariance matrices suggests 
that at least two to three tim.es this minimum algebraic 
requirement of d+1 samples is preferable [7, 11]. Frequent- 
ly the number of available training samples for such complete 
estimation of the individual class covariance matrices is 
inadequate . 

Two common methods of correcting this difficulty involve 
somewhat arbitrary assumptions about the nature of the class 
covariance matrices. One possible assumption is that all 
classes share the same covariance matrix S, v.’hich is then 
obtained by averaging the individual class covariance 
matrices S . : 
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S = 



K. 

1 



K. 

1 



m 

E 

i=l 



K. 



1=1 1 



j=i 






(III-l) 



A second method assumes statistical independence of the 
features and makes all off-diagonal elements zero, regard- 
less of evidence to the contrary. This approach preserves 
additional detail about individual class structures at the 
expense of disregarding highly correlated features. It 
requires only the calculation of a single feature class 
means and variances, as indicated in Equations A-1 and A-2. 
Although such assumptions are almost surely incorrect, they 
often lead to better classifier performance than a maximum 
likelihood estimator [7] . 

Small sample size also leads to complications when one 
attempts to estimate the error rate of the classification 
algorithm finally adopted. If one chooses to estimate the 
classifier's error rate from the assumed parametric model, 
one finds the result is optimistic to the extent that the 
training samples are peculiar and unrepresentative [21] . 

An empirical approach to determining error rate avoids 
this problem by running the classifier on a set of test 
samples. \*71iere one is faced with a small number of available 
signals for which ground truth is accurate, the confidence 
which may be placed in the results of the empirical scheme 
is marginal. For example, if two errors are made in ten 
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test samples, one may predict with 95% confidence that the 
true error rate of the classifier lies somewhere between 
five and fifty-five percent. The figure of ten test sam- 
ples is not unrealistic if one uses different design and 
test samples to avoid the hazards of "testing on the training 
data." 

B. INDETERMINANCE OF CLASSES 

The absence of a predetermined number of classes presents 
additional problems which must be addressed in any practical 
identification scheme. The first problem is that of estab- 
lishing the criteria for excluding a given sample as a member 
of any of the previously determined classes. The second 
problem is one of associating or "clustering" the resultant 
"unclassified" signals into possible new classes. 

Within the formal structure of minimal risk Bayesian 
classification, one may introduce the concept of an addi- 
tional class v;hich corresponds to the heuristic category 
of "don't knov7." More frequently, this category is inter- 
preted as that portion of the feature space for which no 
feature-conditional class probability exceeds a given 
threshold. That is: 

p(C^|x) < X for all i (III-2) 

An alternative method of establishing this threshold 
for multivariate normal distributions is to specify a 
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maximum Mahalanobis distance from the established class 
means . 



as the boundary for v/hich a sample will be included in a 
given class . 

The second problem, that of clustering the "unclassified" 
signals, is less well defined from the standpoint of Bayesian 
analysis. It represents one of a group of problems known 
as unsupervised learning. Jarvis and Patrick [12] present 
several techniques by which such clustering may be performed 
and illustrate the advantages in graphically displaying 
the clustering process. 

One such clustering technique used in several identifi- 
cation systems is illustrated for a two-dimensional case 
in Figure 2. Thresholds based on Mahalanobis distance are 
established for each cluster. If the Mahalanobis distance 
for a given sample, such as exceeds the outer threshold 

distance from the mean of any known cluster, that sample 
is considered to represent a new cluster and possibly a new 
signal originator. For this new cluster, arbitrary dis- 
tance thresholds based on "typical" class covariance are 
assumed. Subsequent samples (S 2 ) lying within the inner 
threshold are used to update the class statistics; those 
samples (S^) lying between the two thresholds are tentatively 
associated with the class, but do not update the class statistics. 
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FIGURE 2. Guard Zone Clustering 
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It is during the clustering procedure that the man- 
machine interface becomes a critical factor in the pattern 
recognition process . An equally important factor in the 
efficiency of the clustering algorithms is a technique for 
feature selection which provides an adequate basis for 
representing new classes. The initial investigator has the 
obligation to use a training set which as nearly as 
possible represents the entire spectrum of possible feature 
parameters . 

C. INACCURATE GROUND TRUTH 

The difficulty of obtaining accurate correspondence 
between the elements of the initial training set and known 
signal originators represents the most persistent and per- 
plexing dilemma in the design of pattern recognition devices. 
Those signals for which there is a high degree of probabili ty 
that the originator can be determined from external infor- 
mation may not adequately represent the entire spectrum 
of signals. Conversely, the originators of a more 
representative group of signals may not be known. 

One may attempt to choose self consistent classes 
through the use of a clustering algorithm and the applica- 
tion of probability estimator to information confirming 
class membership, a technique best described as "learning 
with a probabilistic teacher." This approach, discussed 
more fully by Cooper [4] , is still rather exploratory in 
nature and is implemented at considerable cost in timeliness 
and computational efficiency. 
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since the empirical methods by which any new classifi- 
cation performance is evaluated must use as ground truth 
the identifications of an older system, indications of any 
improved performance in the new system are perforce intui- 
tive . This is the crux of the circular dilemma v;hich arises 
v;hen one is forced to estimate not probability of error, but 
probability of disagreement. In the event of disagreement, 
which system is right? 



20 



IV. FEATURE SELECTION IN WAVEFORM ANALYSIS 



One important question regarding the process of waveform 
identification has not yet been addressed. All of the tech- 
niques described so far have indicated that one commences 
this process with a particular set of feature measurements. 

In the case of a waveform z(t), a sufficiently general 
feature vector might be the sequential sairiples of the 
continuous signal, i.e. {z (t^) , z (t 2 ) • • . , z (t^) } • If these 
samples occur at greater than the Nyquist rate, one may 
reasonably assume that they constitute an adequate descrip- 
tion of the signal. One could then proceed directly to 
the calculation of n-dimensional probability distributions 
for each class, blithely oblivious to the fact that n 
may be on the order of 500 or more. 

In light of the advantages in the choice of a small nuiribe 
of mutually-independent well-separated features for the 
process of pattern identification, it is normally advisable 
to seek some m.apping, F, of the signal space Z into a 
pattern space X. That is, we wish to find a preprocessing 
transformation so that the pattern space X = F(Z) satisfies 
the following objectives: 

(1) low dimensionality 

(2) retention of sufficient information 

(3) enhancement of distance in pattern space as a measure 
of the similarity of physical patterns 
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(4) comparability of features among samples. 

Note that 1 and 2 imply elimination of redundant information. 

In addition to the above objectives one might also wish 
to obtain a pattern space representation adequate for the 
construction of new classes in an unsupervised learning 
mode. To this end, it is convenient that the individual 
components of the pattern space have some natural interpre- 
tation which might provide qualitative information about 
the underlying causes of class difference. 

As a first step toward eliminating redundant informa- 
tion in a pattern space satisfying the above criteria, it 
is frequently desirable to represent the sampled waveform 
as a linear combination of orthonormal functions . Tv;o 
commonly used orthonormal expansions satisfying slightly 
different training set mean square error criteria are dis- 
cussed in Appendix C. The coefficients of these functions 
may then be used as features of the signal. The choice of 
orthonormal functions used in the expansion hov/ever, need 
not be one v;hich necessarily minimizes a form of mean 
square error in the representation of the signal. If one 
is V'/illing to accept a slightly greater amount of error in 
modeling the training set waveforms, then one possible 
technique is to use a set of orthonormal functions for 
which it is computationally efficient, as it is with the 
fast Fourier transform, to obtain coefficients. By trun- 
cating the number of coefficients one m.ay achieve somie 
reduction of dimensionality. 
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The features derived from the orthonormal projections 
of the sampled signal described in Appendix C represent 
linear transformations of the sampled waveform, data. For 
the purposes of signal identification, this representation 
may not be adequate. 

Since one of the goals in this application of pattern 
recognition techniques is to assist a human analyst in 
making an identification, it is particularly desirable that 
the resultant features derived from preprocessing have some 
degree of natural interpretation. The optimal linear trans- 
formations discussed do not alv/ays lend theiriselves to easy 
extension in an unsupervised learning mode. 

Features V7hich are more directly related to an assumed 
model of the signal generating process and v?hich are often 
easily implem.ented as a mieasurement procedure, may in fact 
involve non-linear transformations of the waveform data. In 
selecting features of this type, the prior experience of 
the data analyst and the details of the signal generating 
model serve as a focus for the development of non-linear 
preprocessing. 

Quite commonly, the initial feature space of the signal 
will include components derived from readily interpretable 
data transformations, both linear and non-linear in nature. 
This feature set is normally too large for convenient compu- 
tation and some of the features thus obtained may be depen- 
dent. The techniques of linear dimensionality reduction 
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described in Appendix C may be applied to this set of 
features, particularly v^hen graphical display of the indi- 
vidual signal is required. In order to provide a feature 
set with dimensionality low enough to insure computationally 
convenient classification and general enough to allow the 
establishment of new class clusters, it is both sufficient 
and desirable to reduce dimensioneility while retaining the 
identity of individual features by selecting those n of 
m features v;hich lead to the optimal separation of classes. 
The selection process then requires both the adoption of a 
criterion by vdiich class-separability may be estimated and 
the development of a feature search algorithm. Both of 
these topics are discussed in Appendix C. 



32 



V. PARAtffiTER MAPPING OF RZuSTER SCAN DISPL/vYS 



As a practical application of the signal classification 
techniques previously described, the author conducted a 
series of investigations v/hich used as a data base the raster 
scan displays of bauded, frequency-shift-keyed (FSK) sig- 
nals transmitted by several originators. The signal data 
base, many of the associated processing algorithms and the 
processing equipment itself were made available through the 
cooperation of Electromagnetic Systems Laboratories in 
Sunnyvale California. The data base and supporting facili- 
ties are the results of the company's development of the 
Parameter Encoder in a Navy-sponsored program for signal 
measurement and identification. 

The objective of the author's investigation v/as to iso- 
late from the raster scan data a small set of features which 
might prove useful in automated signal identification. As 
a part of the existing clustering and identification tech- 
nique, the raster scan pattern of a given signal was pre- 
sented to the system opei'ator for visual comparison with 
patterns of previously identified signals. It was claimed 
that such comparison was of value in the final stages of 
the identification process, when the classification of the 
signal in question had been narrowed to a few possibilities. 
Display of the information necessary to produce the required 
number of raster scan patterns, however, consumes a significant 
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amount of computer time. The interests of improved effi- 
ciency, coupled with the author's desire to quantify a 
process open to considerable latitude in operator judgment, 
provided the motivation for the project. 

In obtaining a reasonable initial set of parameters for 
consideration by this thesis, a study of the signal generating 
process and the characteristics of the raster scan display 
led to a slight modification of the normal minimal phase 
raster representation. The resulting "phase-unwrapped" 
raster display v\?as the basis of a series of measurements 
related to clockrate, transient x^h^nomena, and possible 
data- dependent or unintentional external rate modulation. 

A model of the signal process, the raster scan display, and 
the m.ethods used to measure the initial set of features 
are discussed in the follov;ing sections. 

A. SIGNAL GENERATION PROCESS 

Since the raster scan display is based essentially on 
zero crossing information, it represents a non-linear trans- 
formation of the incident signal. It is convenient to refer 
to a model of the signal generating process in an atter:ipt 
to justify use of such a transformation. From this m.odel 
the qualitative effects of the various steps in the genera- 
tion of the raster scan pattern can be estimated. Fig. (3) 
provides an outline of the process which leads to the raster 
scan display. All of the signal data base information avail- 
able for the feature measurements prox^osed by this thesis is 
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contained in the transition time measurements provided by 
the edge detector. 

Of all the qualitative effects upon transition time 
measurements, those associated v^ith the originator's clock 
and bit-stream generator are potentially of greatest value 
to a raster-scan-based identification scheme. Class-to- 
class differences in the basic clock rate are certainly of 
interest, as are clock rate variations due to component 
instability and unintentional external rate modulation. 

If the originator's clock is not continuously running, 
transient phenomena associated with clock turn-on may be 
apparent in the raster display. The binary bit stream 
generator (data modulator) and non-linearities in the F3K 
modulator itself may produce an unequal mark-space duty 
cycle or data dependent mark-space asymmetries which appear 
as bias in the raster display. Such effects should produce 
features of the raster display which serve to characterise 
the originator. To isolate them effectively, how’ever, may 
require the use of non-linear waveform daita transformiations . 

Other parts of the signal generating process may intro- 
duce qualitative effects which tend to obscure those produced 
by the source. For example, noise and signal attenuation 
in the propagation path may introduce FM spikes and anoma- 
lous transition times. Non-linearities and zero-level off- 
set in the FM discriminator may produce additional bias in 
the measured mark and space durations. The waveform, sampling 
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procedure itself introduces an unavoidable time quantization 
error v/hich may or may not be uniformly distributed on a 
short term basis. If this error is uniformly distributed, 
one can expect a transition time error variance v;here 

a = ~ V-T 

12 *■ 

and where At is the sampling period. 

B. RASTER SCAN DISPLAY 

The use of a raster scan display arises rather naturally 
in the attempt to represent the fine details of clock rate 
variation. The axes of the display correspond to coarse and 
fine divisions of time in the following manner. Suppose 
some event, such as a mark-space transition, is assumed to 
occur only at times characterized by some integer multiple 
of a nominal period T. The actual time of occurence, t, may 
then be represented by 

t = nT H- z 

-T/2 < z £ T/2 ; n integer V-2 

By equating the horizontal axis of the raster display to 
the actual time variable, t, and the vertical axis to the 
fine time z, any arbitrary time may be represented as a 
point on the sawtooth waveform indicated in Figure 4b. The 
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FIGURE 4. Relation of Phase Perturbation to Raster Display 
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raster display has the effect of mapping the points t^ 
in a one-dimensional region into the points (z(t^),t^) 
of a two-dimensional region. Since the fine time variable, 
z, represents a lag or lead relative to the nominal period, 
one can relate the points of the raster scan display to 
the samples of some continuous phase perturbation . 

For example, consider the continuous phase function 
illustrated in Figure 4a. The function is comprised of a 
linear part and a perturbation 

<^(t) = 2TT(i t + z(t)) V-3 

where T is the nominal period of the raster scan display 
at times t^. If samples of this function are taken when 

<f>(t) = 2mr n=0,l,2... V-4 

n 



and if 



|z(t) 1 < T/2 for all T V-5 

then the minimal time difference from the expected sample 
time is proportional to the phase perturbation. The raster 
scan points of Figure 4b represent non-uniform samples of 
this continuous perturbation. In this sense one can repre- 
sent the events occuring at times t 2 ...) as samples 
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(z(t^) z(t 2 ) ...) of a continuous phase proportional function 

z(t) where 



z(t„) = 



n 



n 



max 



T 



where 



n 



max 




V-6 



i.e. the greatest integer less than . 

Figures 5a through 5d illustrate the raster patterns 
of events t^ of arbitrary starting phase and periodicity 
different from the nominal period of the raster scan. If 
the nominal raster period coincides with the actual periodic- 
ity, as in Figures 5a and 5b, the slope of the line joining 
points of the raster pattern is zero. Similarly, if the 
periodicity of the series of events differs slightly from 
the nominal raster period by a small amount e, the rastej; 
pattern may be characterized by a line whose slope is given 
by: 



^2 “ ^1 _ e 
^2 “ ^1 T + e 



V-7 



as shown in Figures 5c and 5d. 

In the signal display of the Parameter Encoder both 
mark-space and space-mark transition times are plotted on a 
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raster scan whose nominal period is approximately one-half 
the time for a single mark-space cycle, that is raster 
period is equal to the baud length. The alternate points 
of the raster scan exhibit a bias term if the mark and 
space durations are not equal. If the period of the entire 
mark-space cycle remains constant, lines connecting alter- 
nate points of the raster scan have the same slope, but they 
are vertically offset by the difference in mark and space 
durations. Figure 6 illustrates this condition. 

The effect of the uniform sampling of the waveform and 
the subsequent time quantization of events is apparent when 
the nominal period of the raster scan is on the order of 
ten times tliat of the sampling period. Figures 7 and 8 
illustrate the effect of quantization noise in tv/o raster 
scan patterns of a slowly varying signal, one of V7hich 
uses a nominal period that is is an integer multiple of 
the sampling time. This quantization phenom*ena was evident 
in the raster displays of the signal data base used in this 
thesis since the nominal baud length is on' the order of 
tv7cnty times the sampling period. 

The abrupt transition in Figures 7 and 8 also serves to 
illustrate the "aliasing" effect which may occur in a minimal 
phase representation if the timie difference between the 
actual and the nominal event exceeds one-half the nominal 
period. When this aliasing occurs the points plotted on 
the raster scan are no longer proportional to the assumied 
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FIGURE 6. Raster D5.splay of Bauded Signal Bias 
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FIGURE 7. Quantization and Aliasing Effects in Raster Display 
Nominal Period = 18.6 x Quantization Interval 



Perturbation = A sin 2TTft 
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FIGURE 8. Quantization and Aliasing Effects in Raster Display 
Nominal Period = 20 Quantization Units 



continuous perturbation. In order to maintain proportionality 
it is preferable to represent the time of events by a minimal 
phase different plot. Such a procedure performs a certain 
amount of "phase unwrapping" by assuming that the absolute 
raster time difference between two successive points never 
should exceed T/2 . Figure 9 illustrates the effect of this 
minimal phase difference algorithm. Since each segment of 
the sawtooth pattern raster scan may be extended in either 
forward or backv/ard in time the minimal phase difference 
representation uses this extension to represent those points 
for which the minimal phase difference is less than T/2 . 

In Figure 9, those points occurring in the heavily lined 
region indicate the normal raster scan pattern v/hile a 
continuous line joins the points of the minimal phase 
difference pattern. Since some of the "aliasing" associated 
difficulties in the measurement of features may be eliminated 
by use of a minimal phase difference, this representation 
was included as the initial step in the raeasurement process. 
Figure 26 outlines the flow of the algorithm actually used 
in the measurement program. 

C. RASTER SCAN DESCRIPTIVE PARAMETERS 

The first efforts of the author's investigation of the 
raster scan data were directed towards obtaining a rather 
large set of descriptive parameters, particularly tliose 
which seemed useful in the quantitative representation of 
the effects anticipated from a study of the model signal 
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aliased raster point 

FIGURE 9 Minimal Phase Difference Raster Display 



generating process. Manipulations of the transition time 
data to provide information about the nominal clock rate, 
it's transient variations and external modulations, as well 
as source-associated mark-space bias and data-dependent 
transition time anomalies were the areas of particular 
interest. The parameters considered included representa- 
tions of all or part of the signal's raster phase by poly- 
nomial mean square fits, Laguerre polynomials, and Fourier 
power spectrum terms, as well as measures of m.ean bias, 
intrinsic and total signal variance. The techniques 
employed to calculate these param.eters are described in 
the follov;ing sections. 

]. . Mean Bias and the Effect of Transients 

Although several sources of apparent mark-space bias 
other than those associated with the originator have been 
indicated in the study of the process model, the measure- 
ment of this parameter may provide an adequate estimate of 
the source-associated component. More recent studies of 
new signal data for v^hich the bias effects not associated 
with the originator are believed to be negligible,- indicate 
that suclj source-dependent bias not only exists but can in 
fact be time- varying. 

The simplest estimate of mean bias may be obtained 
by the straightforv.^ard averaging of the raster phase differ- 
ence of alternate points on the raster scan. For the raster 
scan display of the alternate positive and negative transitions. 
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z (t^) , z (t 2 ) / z (t^) . • . / z / the mean bias b is 

given by: 
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V-8 



which is the difference in the mean raster phase of the two 
bias conditions. 

One difficulty which arises from this representation 
of mean bias occurs as a result of clock transient behavior. 

If at any point t^ , raster phase z(t^) differs from the 
raster phase of the starting point by an amount approxi- 
mately equal to one-half the nominal raster period, the 
cyclic representation (minimal phase) of the raster display 
introduces an "aliasing" error into the mean bias estimate 
of Equation V-8. 

To eliminate such errors in the bias estimate, one 
must use the minimal phase difference representation discussed 
previously for each bias condition. 

If the transient behavior of the clock rate is suffi- 
ciently violent, even the minimal phase representation may 
not be adequate. For example, assume that the instantaneous 
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originator clock frequency f, is characterized by turn-on 
from a dead stop: 

f = f [1 - exp (-at)] V-9 

ss 

where a is the reciprocal of turn-on constant, t. The 
instantaneous phase (})(t) in cycles is obtained by integrating 
f from 0 to t; 

f 

= {f^^ t t — — [-1 + exp (-at)]} V-10 

S S 3, 

If transitions occur when 4)(t) = 0,1,2,..., (Figure 10a), 
the raster scan points of Figure 10b are obtained. It should 
be noted that the transition times are non-uniform samples 
of the actual phase. Similar computations may be performed 
for the raster phase of a clock whose transient frequency 
)3ehavior is characterized by a shift from a free running 
frequency, f^, to the steady state frequency, fgg* 

V-11 
V-12 

The above results suggest that the raster-phase 
pattern or a minimal phase difference pattern may be 



f^ = f^ exp (-at) + f^g [1 - exp (-at)] 

= ® *ss 

4>(t) = [exp (-at) -1] + fgg t 
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extremely misleading during the transient period. It is 
evident that the predicted raster phase (in cycles) may 
differ from the actual steady state phase by an amount 
equal to 



f - f 
o ss 



T - T 
ss o 

T T 
o ss 



which may be considerably greater than one cycle. For 

example, if T = .9t and x = 20 T then difference in 
^ o ss ss' 

steady state phase may be more than two complete cycles. 

Using the same time constant x, turn-on from a dead stop 
(f^ = 0) results in a steady state phase difference of 
twenty cycles (from nominal) . The raster scan disp3.ay of 
Figure 11 serves to illustrate tlxe transient turn-on effects 
observed when x == 10 T In such a case the anti-aliasing 
capability of the phase unwrapping process is severely 
strained and it V'^ould be wise to avoid measurement of mean 
bias in this transient region. 

2 . Pattern Variance and Intrinsic Variance 

Pattern variance and intrinsic variance measure 
the extent to v;hich the raster pattern is explained by the 
nominal clock rate and the minimal phase difference represen- 
tation. Large values of total pattern variance indicate 
that the nominal raster chosen does not completely explain 
the pattern. Large values of the intrinsic variance suggest 
that a minimal phase difference representation is not 
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FIGURE 11. Effect cf Transient on Normal Raster Display 



adequate, perhaps as a result of violent transients or 
noise. 



These two measures are averaged over points of both 
bias conditions and z^^^ corresponding to the alternate 

positive and negative going transitions. The total variance 
v^ is obtained from the follov;ing calculation: 
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Intrinsic variance (v^) is estimated by; 
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where refers to the jth point of the set of either 

positive or negative going transitions. 

These two terms may be particularly useful as a 
measure of signal degradation. They may be used as a measure 
of non-parametric correlation since the fraction of total 
variance "explained" by the intrinsic variance is given by; 



A(z) 




V-16 
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One estimate of non-parametric correlation is: 



B(z) = Va(z) . V-17 

Small values of B(z) suggest that the apparent behavior of 
the raster pattern is not well characterized by a relatively 
smooth function of time. 

2 • Polynomial Mean Square Fits 

In the hope of quantifying raster pattern behavior 
v;hich system operators had described as "curve up," "curve 
down," flat," or complex," a fourth order polynomial mean 
square fit to the minimal phase difference raster data over 
the entire signal period v/as attempted for each bias condi- 
tion. The fourth order fit was initially believed to be 
capable of representing both transient and steady state 
phenomena. Unfortunately, even the minimal phase difference 
representation of the data could not eliminate aliasing 
problems in those displays where apparently large initial 
transients occurred. As an alternative, a second order fit 
over that data lying outside the region of worst transient 
behavior was eventually incorporated into the measurement 
program. 

The coefficients of the second order fit also serve 
the purpose of reducing the sensitivity of the Fourier trans- 
form data to phase offset and endpoint discontinuities which 
result from minor adjustments in the nominal raster period 
of the display. 
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As an example, consider the periodic Fourier series 
transform of a pattern characterized by a line with non-zero 
slope. As a result of endpoint discontinuity in a periodic 
representation, the odd Fourier series components exhibit 
amplitudes proportional to the magnitude of the discontinuity 
and decreasing in frequency as 

It has previously been established that transient 
phenomena, zero point offset, and differences between the 
nominal and actual steady state period can produce patterns 
V7ith both offset and non-zero • slopes . If the raster pattern 
shows significant curvature, it is difficult to justify any 
particular choice for the nominal period. Any choice of 
nominal period creates its ov;n contributions to the Fourier 
series components. Since the nominal raster period of the 
display is generally the result of an automatic measurement 
procedure, and since it represents one of the signal features 
currently used by the system's identification algorithm, the 
interests of comparability were best served by partially 
correcting for these discontinuities in terms of the m^ean 
square fit coefficients. Only the difference betv?een the 
raster pattern and the polynomial mean square fit to each 
bias condition was used as input to the fast Fourier transform 
calculation. 

4 . Fourier Power Spectrum 

In order to observe the effects of possible uninten- 
tional external modulation and periodic variation of the clock 
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rate, the use of Fourier power spectrum components was pro- 
posed. The calculation of these terms presented an intriguing 
problem, since the minimal phase raster data represent non- 
uniform time samples of a presumably continuous signal. The 
non-uniformity of sampling tim.e arises from both the clock 
rate variation and the fact that the transition times them- 
selves and data modulated. The presence of one transition 
per baud is not the case; transitions occur at various 
integer numbers of originator clock cycles. 



ous if one intends to obtain only a fev; Fourier power series 
components. The techniques for estimating amplitude phase 
and frequency of the principal Fourier components has been a 
topic in recent literature [ 22 ] . One can ea.sily construct 
Fourier series of terms which are orthogonal on the interval 
[0,L] by considering only the frequency terms which are 
periodic in L. 



The problem of non-uniform sampling is not too scri- 



z(t) = Z c^e 



oo 




V-18 



These complex components 
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/ z (t) e 




dt 
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0 



may be approximated step-wise 
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7 , z (t^ ) c 
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V-2 0 



c 



n 



or by numerical integration involving higher order interpo- 
lation. Although numerical integration of the above type is 
possible, all frequency components up to the nominal clock 
frequency were potentially of interest as features . An 
alternate form of calculating the magnitudes of the pov;er 
spectral components, using the fast Fourier transform V7as 
adopted for computer implementation. 



used about 75 non-unif ormly spaced raster data points obtained 
outside the transient region. For each condition, 128 
sample points uniformly spaced over the same time interval 
were chosen. The raster phase at these points v/as obtained 
by an interpolation procedure v/hich performed a cubic poly- 
nomial fit to four minimal phase difference raster points, 
two on either side of the interpolation point. Since the 
fast Fourier transform routine used accepts complex values 
as input, the interpolated functions for positive going and 
negative going transitions were stored as the respective 
real and imaginary parts of this input. The complex input 
function may be expressed as: 



For each bias condition, this procedure typically 



y = f + ig 



V-21 
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where the nth Fourier component of the real ftmctions f and 
g are: 



F^(f) = a^ + i 
n 11 

Fn(g) = a2 + i b2 V-22 

The linearity of the Fourier transform and the requirement 
that for a real function 

F ^(f) - F*(f) V-23 

-n n 

implies that the nth Fourier transform components of both 
functions f and g may be recovered from (F^(y) and F_j^(y)) / 
the real and imaginary Fourier components of the complex 
input [10] . Additionally it can be shovm that: 

[Re[F^(y)]]^ + [Im[F^^(y) ] ]^ + Re (F (y) ) ] ^ + [Im [F_^ (y ) ] ] ^ 

= 2(a^^ + b^^ + a2^ + b2^) V-24 

The Fourier transform information used as features in this 
thesis was the squared magnitude of the frequency components 
summed over both bias conditions in the manner shov;n above. 

The Figures 12-22 show the behavior of the time 
quantization interval, the minimal phase difference interpola- 
tion, the second order polynomial fit, and the Fourier transform 



59 



features on a signal whose possible transition times t^ 
are given by 

t^ = {nT + Asin 2iTfnT + B exp(-anT)) V-25 

The data points were selected from these possible transition 
times by a simulated data modulation which repeated itself 
about times in the duration of the display. In the ras- 
ter patterns the values of the second order polynomial mean 
square fit and the interpolated raster phase for one bias 
condition at the 128 points used for interpolation are shovm 
as the characters (+) and {*) respectively. 

The effect of the quantization interval is apparent 
in all raster displays, but it is encouraging to note that 
a 25 Hz frequency component of amplitude equal to one-half 
the quantization interval (Figure 20) may be detected even 
in the presence of a transient during the initial half second 
of the signal sample. This detection ability, hov/ever, is 
strongly dependent on the apparently uniform distribution of 
transition times within the quantization interval. For this 
deterministic case, as the raster period approaches an integer 
multiple of the quantization interval, the ability of the 
measurement technique to isolate a signal of small amplitude 
as a single maxima deteriorates. 

5 . Laguerre Polynomial Coefficients 

In the attempt to characterize the apparent transient 
behavior of some of the raster displays, a set of Laguerre 
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FIGURE 12. Polynomial Mean Square Fit and Interpolation, Signal 
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FIGURE 13. Fourier Magnitude Components - Signal 



T = 18.6 \inits 
Perturbation = A sin 2irft 
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FIGURE 14. Polynomial Mean Square Fit and Interpolation, Signal 
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FIGURE 15. Fourier Magnitude Components - Signal 



Perturbation = A sin 2iTft 
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FIGURE 16. Polynomial Mean Square Fit and Interpolation - Signal 
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FIGURE 17. Fourier Magnitude Components, Signal 
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FIGURE 18. Polyncmial Mean Square Fit and Interpolation, Signal 
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T = 18.6 units | 

Perturbation = A sin (2iTft) • + B exp{~at) 
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FIGURE 20. Raster Display, Signal 
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FIGURE 2i. Polynomial Mean Square Fit and Interpolation, Signal 
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FIGURE 22. Fourier Magnitude Components, Signal 



polynomial coefficients were calculated over all the points 
z'(tj^) of the initial portion of the signal. The coefficients 
of the first five Laguerre polynomials were obtained by 
stepwise numerical integration over the first 800 quantization 
intervals . 



Lf. = Z z'(t.) exp [-pt . ] ( At . ) 

^ t.<800 

1 — 



= Z z' (t . ) [2pt .-I] exp [-pt . ] (At. ) 



= Z z' (t. ) [2 (pt . ) ^-4pt.+l] exp [-pt . ] (At. ) 

tj ^<_800 111 11 



Lo = Z z ' (t . ) [4(pt . ) ^" [6pt . ] ^+6pt . -1] exp [“Pt. ] (At . ) 

t ^£800 1 .5 1 1 1 11 



L. = Z z' (t . ) [-|(pt . ) '^-~(pt. ) ^+12pt .-8pt .+1] exp [-pt . ] (At . ) 
t.<800 1 . 51.31 1 - 11 

^ V-2 6 
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Ati = t.^3_ - t. 



During the transient period, the aliasing effect, 
which occurs if the nominal raster period does not match 
the actual period, is especially pronounced. The data may 
very well be undersampled in this region, thus no attempt was 
made to compute these coefficients independently for each 
bias condition and the phase unwrapping procedure was applied 
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to all data points. This process, hov;ever, operated on the 
premise that the phase difference from transition to transi- 
tion is not more than one-half the nom.inal raster phase. 

Even if both positive and negative transitions are used for 
phase unwrapping, this premise may not be valid. 

If the initial transition times predicted by the 
steady state period and knowledge of the signal’s underlying 
data modulation can be determined, a more accurate represen- 
tation of transient phase might be obtained. Unfortunately, 
knowledge of this data modulation v;as not available for the 
purpose of this thesis. 

D. COMPUTER FILES AND PROGRAMS 

The programs used to produce the feature measurements 
described in the previous section were designed to operate 
on a pre-production model of the Parameter Encoder which 
incorporated in its hardv^are a Hewlett-Packard 2100A computer 
with tv;in disc drives and moving head disc operating system. 
The normal graphics display and list device was a Tektronix 
storage scope. 

The programs discussed in this section, v^hile written 
in HP Fortran, make use of the Fortran-callable executive 
routines available under the disc operating system. These 
routines permit program overlays and disc file input/output. 
Similarly, a library of Fortran-callable utility routines 
developed at Electromagnetics Systems Laboratory to control 
graphics in put and output were used to create and erase the 
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scope displays, temporarily halt computations and provide 
keyboard control of the measurement process. 

The signal data base used as input to the measurement 
process was stored in the disc file, PART. This file con- 
tained a signal index number, the automatically measured 
clock rate, the time difference between successive signal 
transitions and the total number of transition times 
measured (see Table 1) . Since this data format differed 
slightly from that normally accepted by the standard raster 
display routine RASF, an interface program FACE was v^ritten 
to convert the data. 

Following the initial display of the raster signal, 
process control could be transferred to a program. FOURD 
v/hose function v;as to call the interpolation, measurement, 
and display subroutines and to store the measured parameters 
in the disc file PAR. The overall process flow of FACE, 

RASF and FOURD is shown in Figures 23 and 24 . 

1 . Measurem.ent Process 

The first subroutine called by FOURD is RESDU (see 
Figures 25 and 26 for overall process flow ) . This subrou- 
tine calculates a minimal phase difference representation 
for the points of both bias conditions outside the transient 
region. The points thus represented are used by the subrou- 
tine CURV to calculate the coefficients of the polynomial 
mean square fit. They are also used to calculate the signal's 
mean bias and its total and intrinsic variance. A cubic 
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TABLE I 



Disc File PARF 



Size: 1542 sectors - 128 Integer words per sector 



Sectors 0-5 Directory - 768 words 



Word Use Comment 

1-768 Signal number Order in directory 

corresponds to order 
of signal files 



Sectors 6-1542 Signal Files 2 Sectors per signal 



V7ord 


Use 


Comment 


1-2 


Nominal period 


Floating point 


3 


Number of 
measurements 




4 - 256 


Transition 

time 


Time difference between 
successive transitions 
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FIGURE 23. Process Flow, Interface and Raster Display 
PrograiTts FACE and RASF 
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FIGURE 




24. Process Flow Measurement and Display 
Program FOURD 
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FIGURE 25, Flow Diagram Interpolation and Measurement 
Subroutine RESDU 
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FIGURE 26. Minimal Phase Difference Algorithm 
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interpolation scheme CURF then provides the uniform samples 
used by the fast Fourier transform subroutine. RESDU also 
calls the subroutine LAGER, which performs a stepwise numeri- 
cal integration to estimate the coefficients of the first 
five Laguerre polynomials . 

The two other subroutines called by FOURD are F0UR2 , 
a standard fast Fourier transform, and FORD, a routine v;hich 
calculates the average square magnitude of the Fourier com- 
ponents and displays the results on the Parameter Encoder's 
storage scope. The disc file, PAR, v;hich contains all the 
above feature measurements, is organized as shown in Table 1.1 . 

2 . Statistics and Formatting Routines 

The analysis of the features measured by the routines 
indicated in FOURD makes extensive use of identification and 
analysis software previously developed for the Parameter En- 
coder. These procedures use a yet another disc file structure, 
MASTF (Table IIH , as their input data base. Additional refor- 
matting and selection of the parameters was necessary to allow 
the use of this file since it accepts a maximum of 50 
parameters . 

The program CSTAT provided class statistics for the 
measurements stored in MASTF. In addition to reformatting 
individual signal parameters, and storing them in the disc 
file MASTF, the subroutine TRNSG performed the important 
function of interpolating the existent Fourier magnitudes 
to estimate the square magnitude components at 2 Hz incre- 
ments. This was necessary since the resolution of each signal 
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TABLE II 



Disc File PAR 

Size: 200 Sectors (128 Integer VJords per Sector) 

Sector 0 Directory 

Word Use Comment 

1-128 Signal Number Order in directory corres- 

ponds to order of files 

Sectors 1-199 Signal Parameter Files 2 Sectors per File 



Word 


Use 


Comment 




1-2 


Nominal Period 


Floating Point 


SIG (1) 


3-4 

5-8 


Signal Resolution 
(Not Used) 


Floating Point 


SIG(2) 


9-18 


Polynomial Coefficients 
Mean Square Fit 


Bias Condition 


I SIG(5-9) 


19-28 


Polynomial Coefficients 
Mean Square Fit 


Bias Condition 


2 SIG (10-14) 


29-38 


Polynomial Coefficients 
Laguerre 


Floating Point 


SIG(15-19) 


39-40 

41-70 


Mean Bias Estimate 
(Not Used) 


Floating Point 


SIG (20) 


71-72 


Variance 


Floating Point 


SIG (35) 


73-74 

75-78 


Intrinsic Variance 
(Not Used) 


Floating Point 


SIG(36) 


79-208 

209-256 


Fourier Magnitude 
Components 

(Not Used) 


65 Terms SIG(4C 


1-104) 
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TABLE III 



Size: 1524 

Sector 0 

Word 

1 

2 

3 

5 

6 

14 

15-64 



Signal File 

Word 

1 

2-3 

4 

5-104 

105-154 

159-384 

Class File 

Word 

1-2 

3 

4 

5-104 

105-204 

205-254 

255-304 

305-384 



Disc File HASTF 

Sectors 128 Integer V7ords per Sector 
File Information 
Use 

Number of Sectors in MASTF 
Number of Words in Directory 



Number of Words in File 
Scale Code (-1 No scaling) 
Variance Breakpoint 
Number of Parameters 
Parameter Numbers 

Use 

Signal Number 
Identification Nui'ober 
(Not Used) 

Parameters 
Parameter Weight 
(Not Used) 

Use 

Identification 

Number of Entries 

Number of Parameters 

Class Means 

Class Variance 

Number of Parameters Used 
for Classification 

Parameter Numbiers 

(Not Used) 



Comment 



50 Floating Point Words 
1 If used, 0 If Not Usee 



Comment 



Floating Point 
Floating Point 
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was variable. Ground truth classification of each signal 
had been previously determined through the use of parameters 
not included in this thesis, and was stored in a small Disc 
File, CSIGF. 

E. AIJALYSIS OF PARAT4ETERS 

From the basic group of features obtained by the measure- 
ment process, a group of 50 parameters v;as selected for further 
analysis. The analytic techniques used were those readily 
available to the user of Parameter Encoder identification 
software. Several iterations of measurement, observation, 
and feature analysis were conducted before arriving at the 
parameters indicated in Table IV. Earlier measurements had 
involved the use of fourth order polynomial mean square fits 
to the entire raster display, and the use of Fourier magni- 
tude coefficients of raster pattern input uncorrected by 
the polynomial fit. 

These earlier techniques appeared to provide reasonable 
measurements for approximately 60 percent of the signals 
considered, but failed catastrophically for the remaining 
signals. The major cause of failure appeared to be the ina- 
bility of the minimal phase difference representation to 
track the raster pattern through initial transients and 
anomalous transition times. The final measurement process 
discussed in the previous section was applicable to about 
80 percent of the signals attempted and resulted in measure- 
ments for the final training set of 74 signals distributed 
among nine classes as shown in Table V. 
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902 

903 

904 

905 

906 

907 

908 

909 

910 

911 

912 

913 

914 

915 

916 

917 

918 

919 

• 

950 



TABLE IV 



Parameter Identification 
Identification 
Nominal Period 
Length of Signal Resolution 



0 


Order 


Polynomial 


Coefficient 


Bias 


Cond . 


1st 


Order 


Polynomial 


Coefficient 


Bias 


Cond . 


2nd 


Order 


Polynomial 


Coefficient 


Bias 


Cond . 


0 


Order 


Polynomial 


Coefficient 


Bias 


Cond . 


1st 


Order 


Polynomial 


Coefficient 


Bias 


Cond . 


2nd 


Order 


Polynomial 


Coefficient 


Bias 


Cond . 



1st Laguerre Coefficient 
2nd Laguerre Coefficient 
3rd Laguerre Coefficient 
4th Laguerre Coefficient 
5th Laguerre Coefficient 
Mean Bias 
Total Variance 
Intrinsic Variance 

Fraction of Total Variance Explained by Intrinsic 
Fourier Power Component 0 Hz 
Fourier Power Component 2 Hz 



Fourier Power Component 64 Hz 
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TABLE V 



Training Set Signal Distribution 



Class 


1 


(Not Used 


in 


Distance 


Ratio) 




1447 




1300 




875 


12 6 


171 


172 




775 




264 






Class 


2 












All 




758 




676 






Class 


3 












770 




183 




192 


1483 


14 8 8 


693 




1354 




1367 


1498 


818 


306 




307 




319 


1525 


938 


Class 


4 












1441 




1444 




483 


1309 


771 


1327 




1337 




1338 


742 


315 


1399 




925 










Class 


5 








• 




127 




155 




766 


826 




Class 


6 












489 




235 




236 


787 


1543 


Class 


7 












1286 




1289 




672 


475 


674 


675 




120 




1304 


1305 


141 


Class 


8 


(Not Used 


in 


Distance 


Ratio) 




888 




239 




1487 


698 


915 


446 




970 










Class 


9 












121 




754 




145 


149 


196 


189 




201 




202 


203 


204 
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The ground truth identification of these signals was 
based on the agreement of two external identification 

s 

schemes. One of these was the semi-automated process 
developed by Electromagnetic Systems Laboratory v^hich uses 
additional signal parameters not associated with the raster 
display. The other process uses all-source information in 
arriving at signal identifications. 

The analysis began v;ith the measurement of single fea- 
ture separating capability for each of the 50 selected 
parameters. The eleven features showing the greatest capa- 
bility were then examined for redundant features by calcula- 
tion of class and global correlation matrices. Finally, the 
error rate of a minimum distance classifier v;as used as the 
criteria for a feature search procedure v;hich combined single 
feature ranking and search-vfithout-replacem.ent techniques. 
Since the amount of training data v/as limited, even this 
iterative procedure may place too great an emphasis on the 
training set values . 

1 . Single Feature Separability 

The single feature separability measure available 
in the Parameter Encoder softv;are is a modification of the 
distance ratio techniques discussed in Appendix C.2. The 
program FINTR calculates the average square distance between 
points of m different classes. 




2 

m (m-1) 



m-1 m 

Z E 
i=l j= i 




K. K . 
1 3 

Z Z 
k=l 1=1 



(X 



(i) 

1 




) 



2 



V-27 
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The m class average within-class-variance , 




V-2 8 



and uses the ratio 



D. 



2 



F 



B 



V-29 



D, 



w 



as a measure of separability. Since the array storage 
available in this program limits the number of input samples 
to 60, the signals of classes 1 and 8 were not used in the 
calculation of the feature distance ratios. The distance 
ratios obtained using all signals of the other 7 classes 
are shown in Table VI. The distance ratios of the frequency 
components are also plotted in Figure 27. 

By way of comparison Figure 28 shows the distance ratio 
as a function of frequency for data obtained in earlier 
investigations using a smaller number of signals and classes. 
The V7ell defined maxima in the latter case can be attributed 
to the fact that the signals and classes user v;ere particu- 
larly noise and transient free. In these cases one was able 
to use the entire raster display uncorrected by a polynomial 
fit to estimate the Fourier components. 

If one assumes that the parameters are statistically inde- 
pendent, then the ranking of features according to distance 
ratio provides a guide for feature selection. Eliminating 
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§ 

901 

902 

903 

904 

905 

906 

907 

908 

909 

910 

911 

912 

913 

914 

915 

916 

917 

918 

919 

920 



TABLE VI 



ingle Distance 


Ratios 


for Parameters 


; 901 


- 950 


Ratio 


# 


Ratio 


# 


Ratio 


10.59 


921 


1.32 


941 


1.14 


11.68 


922 


1.14 


942 


1.08 


1.87 


923 


1.35 


943 


.94 


1.07 


924 


1.25 


944 


.93 


1.07 


925 


1.45 


945 


.93 


1.17 


926 


1.81 


946 


.84 


1.00 


927 


2.10 


947 


.95 


1.03 


928 


1.51 


948 


.92 


1.15 


929 


1.09 


949 


.96 


1.03 


930 


1.44 


950 


.97 


1.01 


931 


1.68 






.99 


932 


1.06 






1.08 


933 


1.27 






.98 


934 


1.11 






.99 


935 


1.28 






1.20 


936 


1.44 






1.14 


937 


1.35 






1.68 


938 


1.26 






1.22 


939 


1.02 






1.19 


940 


1.01 
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FIGURE 
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27. Feature Separatidn Quality vs. Frequency, 
59 signals, 7 classes 
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FIGURE 28. Feature Separation Quality vs. Frequency, 
30 signals, 4 classes, no polynomial fit 
corrections 
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a 






those features v/hose ratio was less than one and ranking 
the best eleven of the remaining resulted in the feature set 
shown in rank order in Table VII. 

2 . Multiple Feature Analysis 

To test the assumption that the parameters selected 
by single feature ranking were independent a utility routine 
PCORL was used to calculate their class and global correla- 
tion coefficients. Tv/o dimensional scattergrams for selected 
pairs of highly ranked features were used for visual confir- 
mation of feature correlation and clustering. In the global 
case, two of the parameters exhibited correlation coefficients 
in excess of 0.9. These parameters, whose scattergram is 
shown in Figures 29 and 30, are the nominal raster period and 
the signal resolution . Correlation between these two 
parameters might be expected to be high if the sample signal 
length corresponds to the same number of bauds for all signals. 
Both the scattergram and the class correlation coefficients 
indicated that the period and signal duration remain highly 
correlated on a class basis. Other scattergrams such as that 
shown in Figures 31 and 32, for the next most effective single 
feature and the nominal raster period, gave little evidence 
of correlation. The results of these calculations and obser- 
vations for other feature pairs were sufficiently encouraging 
that single-feature distance ratios were used to choose the 
next four features to be discarded from the set of eleven. 



\ 



91 



£ 

902 

901 

927 

903 

926 

931 

918 

928 

925 

936 

930 



TABLE VII 



Ratio 

11.68 

10.59 

2.10 

1.83 

1.81 

1.68 

1.68 

1.51 

1.45 

1.44 

1.44 



Ranked Features 

Identification 
Signal Resolution 
Period 

18 Hz Fourier Component 

0 Order Poly Mean Square Fit, 
Bias 1 



16 


Hz 


Fourier 


Component 


26 


Hz 


Fourier 


Component 


0 


Hz 


Fourier 


Component 


20 


Hz 


Fourier 


Component 


14 


Hz 


Fourier 


Component 


36 


Hz 


Fourier 


Component 


24 


Hz 


Fourier 


Component 
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RESOLUTION 
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NOMINAL PERIOD 

FIGURE 29. Scattergram of Class Membership, Nominal Period vs. Resolution 
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FIGURE 30. Expanded View of Fig 




I 




expanded in Figure 32 




* * % * 
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NOMINAL PERIOD 

FIGURE 31. Scattergram 18 Kz Component vs. Nominal Period 
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FIGURE 32. Expanded View of Figure 31 



96 



The six features finally chosen were selected by 
repeated application of a minimum-normalized distance 
classification program CNFUS . Sixty-two signals from all 
classes except class four were classified by the program 
according to their distance from the class means. This dis- 
tance d^ of a particular set of n signal features, (x^...x^) 
from the class mean ) was given by 



d/ 

1 



1 

n 



n 

Z 

j=l 









(a f ^ 
. D 



V-30 



The signals used for feature selection, their distance 
to the nearest four class means, and the confusion matrix 
resulting from the choice of the nearest class are shov;n in 
Tables VIII and IX. This choice of six features leads to 
approximately 56 percent agreement betv/een ground truth data 
and the first choice of the classifier v;hen the entire train- 
ing set, including class four, is used. This, of course, 
represents an optimistic estimate of the probability of 
correct classification. Further testing of the classifier 
over an independent set of signal data is advisable. That 
the addition of features other than the basic clock rate does 
improve classifier performance is evident in the improvement 
of the 34 percent classification probability obtained using 
only this feature as the basis for classification. The 
resultant single feature confusion matrix is shown in Table 
X. 
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TABLE VIII 



Classifier Selections of 4 Classes Nearest 
to Signals of the Training Set 

Parameters Used 

902 903 927 926 931 918 

Ground 

Signal Truth Classes and Distance to Mean 

Class 



427 


2 


2 


1.16 


1 


1.25 


4 


2.52 


5 


4. 


758 


2 


2 


1.15 


1 


2.80 










676 


2 


1 


.70 


2 


1.12 


4 


1.96 






770 


3 


3 


.57 


8 


.73 


4 


.79 


6 


1. 


183 


3 


4 


.51 


1 


.88 


8 


.94 


3 


1. 


192 


3 


1 


. 55 


3 


.84 


4 


.91 


8 


1. 


1483 


3 


3 


.65 


4 


. 66 


1 


. 84 


8 


« 


1488 


3 


8 


.74 


4 


.90 


3 


.99 


1 


1. 


693 


3 


4 


1.07 


6 


1.23 


8 


1.44 


1 


1. 


1354 


3 


4 


.72 


8 


.77 


3 


.78 


9 


1. 


1367 


3 


3 


1.01 


4 


1.12 


1 


1.15 


8 


1. 


1498 


3 


4 


1.24 


3 


1.25 


6 


1.46 


8 


1. 


818 


3 


3 


.83 


4 


.95 


8 


1.23 


1 


1. 


306 


3 


4 


.40 


3 


.62 


8 


.64 


6 


. 


307 


3 


3 


.76 


4 


.86 


9 


.98 


8 


1. 


319 


3 


3 


.96 


4 


1.14 


8 


1.32 


1 


1. 


1525 


3 


8 


.57 


1 


.68 


4 


.81 


3 


1. 


930 


3 


3 


.54 


8 


.56, 


4 


.84 


1 


1. 


127 


5 


5 


.93 


1 


1.82 


4 


1.88 


8 


2. 


155 


5 


5 


.83 


4 


1.40 


7 


1.51 


6 


1. 


766 


5 


5 


.91 


6 


2.41 


4 


2.53 


8 


2. 


826 


5 


8 


.42 


1 


.57 


4 


.78 


9 


. 


489 


6 


1 


.84 


4 


.89 


6 


.92 


8 


1. 


235 


6 


6 


.88 


4 


.99 


8 


1.61 


3 


1. 


236 


6 


6 


.96 


4 


1.67 


1 


2.23 


5 


2. 



43 

36 

07 

08 

85 

23 

52 

08 

16 

49 

28 

82 

13 

78 

08 

05 

33 

71 

57 

89 

17 

65 

29 
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i 




787 

1543 

1286 

1289 

672 

475 

674 

675 

120 

1384 

1305 

141 

121 

754 

145 

149 

196 

189 

201 

202 

203 

204 

1447 

1300 

875 

126 

171 

172 

775 

264 

888 

239 

1487 

698 

915 

446 

970 



6 


6 


1.21 


1 


1.33 


4 


1.56 






6 


8 


.70 


6 


.73 


4 


.85 


1 


.89 


7 


7 


. 80 


4 


1.03 


1 


1.20 


6 


1.55 


7 


7 


.88 


4 


1.24 


1 


1.41 


5 


1.59 


7 


7 


.88 


4 


.94 


1 


1.00 


9 


1.30 


7 


7 


.61 


4 


1.38 


1 


1.59 


6 


1.76 


7 


7 


1.25 


4 


1.34 


6 


1.51 


9 


1.62 


7 


7 


.95 


1 


.96 


4 


1.57 


5 


1.05 


7 


7 


1.37 


4 


1.59 


1 


1.91 


6 


2.28 


7 


-L 


.47 


4 


.79 


7 


.83 


6 


1.13 


7 


1 


.51 


7 


.65 


4 


1.23 


6 


1.92 


7 


1 


.40 


7 


.99 


4 


1.14 


5 


1.53 


9 


4 


.73 


9 


.89 


7 


1.05 


8 


1.14 


9 


1 


.50 


4 


.53 


7 


.72 


5 


.94 


9 


1 


.65 


7 


.82 


4 


.93 


9 


1.02 


9 


3 


.75 


4 


. 85 


8 


.87 


9 


1.02 


9 


9 


.32 


3 


.59 


4 


.71 


1 


.74 


9 


9 


.63 


4 


. 87 


3 


1.08 


7 


1.09 


9 


8 


. 80 


1 


.80 


4 


. 83 


9 


1.05 


9 


4 


. 56 


3 


.58 


8 


.58 


6 


.83 


9 


9 


1.08 


4 


1.26 


1 


1.41 


3 


1.43 


9 


1 


.93 


9 


1.00 


4 


1.25 


3 


1.32 


1 


1 


1.52 


4 


5.19 


2 


5.24 






1 


1 


1.52 














1 


1 


.73 


4 


2.79 


2 


3.41 


8 


4.68 


1 


1 


.40 


4 


3.75 










1 


1 


.60 


4 


3.22 


2 


4.19 






1 


1 


.79 


4 


3.29 


2 


3.41 


8 


5.64 


1 


1 


.87 


4 


3. 09 


2 


4.12 






1 


1 


.62 


4 


2.56 


2 


2.67 


8 


4.37 


8 


8 


. 61 


1 


.83 


4 


1.30 


6 


1.79 


8 


8 


.55 


3 


.69 


1 


.84 


4 


.89 


8 


8 


.76 


1 


1.16 


4 


1.24 


6 


1.58 


8 


6 


.87 


4 


1.05 


8 


1.05 


1 


1.18 


8 


8 


.90 


5 


.92 


1 


1.24 


4 


1.34 


8 


1 


.66 


8 


1.04 


4 


1.12 


5 


1.29 


8 


4 


1.00 


3 


1.02 


8 


1.32 


1 


1.34 
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TABLE IX 



CONFUSION MATRIX FOR CLASSIFIER USING 
PARAMETERS 902, 903, 927, 926, 931, 918 





1 


2 


3 


CLASSIFIER CHOICE 
4 5 6 


7 


8 


1 


8 


- 


- 


- 


- 


- 


- 


2 


1 


2 


- 


- 


- 


- 


- 


3 


1 


- 


7 


5 


- 


- 


2 


4 


- 


- 


- 


- ■ 


- 


- 


- 


5 










3 






1 


Ul 

(/} 

nj 

•-< 6 

o ° 


1 


— 


— 


— — 


3 


— 


1 


cH 

1 ^ 


3 


- 


- 


- 


- 


7 


- 


8 


1 


- 


- 


1 


1 


- 


4 


9 


3 


- 


1 


2 


- 


- 


1 


# of 


10 


0 


1 


8 0 


1 


0 


5 


e>osess 
Type II 
Error 


GOOD 


37 




MISSED 25 









59% correct classification over training set. 
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# misse 



IVPe 

Error 



0 

1 

8 

0 

1 

2 

3 

3 



3 7 



0 



100 



H a 



TABLE X 



CONFUSION MATRIX FOR CLASSIFIER USING CLOCK RATE ONLY 



1 

1 8 



CLASSIFIER OIOICE 



# missed 
IVpe I 
Error 



m 

<0 

O 6 

7 

8 



# of 
excess 
IVpe II 
Error 



1 

3 



1 

3 

2 



2 

11 



2 

3 



1 

1 

2 



3 

4 



1 

6 



1 10 

d 

3 

4 

4 10 



1 

5 



# GOOD 21 



if MISSED 41 



33.87% correct classification over training set. 



101 



The sequence by which the final six features were 
obtained involved first the calculation of the classifier 
confusion matrix as each of the four lowest ranked features 
were removed from the set of eleven. Their removal improved 
the classifier's performance from 50 percent to 58 percent. 

From this set of seven features, the effect of removing one 
parameter and using the remaining six yielded the following 
classifier performance: 

Parameter Removed Probability of Correct Classif icaticr. 



None 


58 


% 


918 


25 


o, 

'o 


931 


53 


% 


926 


50 


% 


903 


51 


% 


927 


53 


% 


901 


59 


Q. 

“O 


902 


56 


% 



The above suggests that the best single feature to be removed 
at this stage is 901, the nominal raster period. That this 
remaining group of six features is not necessarily the 
optimal choice, can be shown by observing that a choice of 
five features, including parameter 901 leads to a 61 percent 
probability of correct identification shown in Table XI. 
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Actual Class 



TABLE XI 



CONFUSION MATRIX FOR CLASSIFIER 
PAR7\METERS USED 901, 903, 926, 931, 918 



12 3 

1 8 - - 

2 1 2 - 

3 1-6 

4 - _ _ 

5 - - - 

6 1 - - 

7 2 - - 

8 1 -^- 

9 2 - - 

# Excess 800 

Type II 

Error 



CLASSIFIER CHOICE 
4 5 6 7 

4 _ _ „ 

- 3 - - 

^ _ 4 _ 

1 - - 7 

1 1 - - 

3 - - - 

9 10 0 



# missed 

8 9 Type I 

Error 

- - 0 

- - 1 

4-9 
- - 0 

1-1 

1 

- - 3 

4-3 

14 6 

6 0 0 



# GOOD 38 # MESSED 24 



61.29% CORRECT CLASSIFICATION OVER TRAINING SET 
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VI., SUMJyiARY OF CONCLUSIONS AND RECOMMENDATIONS 



It has not been the intent of this thesis to suggest 
feature measurements which can replace those currently used 
in the Parameter Encoder's source identification scheme. The 
fact that the measurement process discussed in this thesis 
is currently xmable to provide meaningful measurements for 
about 20% of the signal data attempted is but one reason for 
not adopting the indicated Fourier components and polynomial 
coefficients as standard features to be measured. 

These features however, constitute a hitherto unexploited 
source of raster pattern information. The use of four of 
these features in conjunction with the nominal raster period 
or signal clock period provides approximately 100% improve- 
ment in classifier capability over the separation performance 
of the nominal raster period alone. In view of the amount of 
overlap in the feature distributions for several classes, 
reflected in both the scattergrams and the distances from 
the individual signals to the nearest class means, it is not 
advisable to use more than four or five such features. Neither 
is it advisable to draw too firm a conclusion as to the opti- 
mal choice of these new features. The 62 test and training 
signals do not by any stretch of the imagination constitute 
a completely representative statistical base. 

It is rather surprising that those terms intended to 
characterize the mean bias, variance measures and transient 
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phenomena show as little single-feature separating capability 
as they do, since it was precisely these phenomena which have 
been used by system operators as visual aids to signal iden- 
tification. While one might reasonably associate the zero 
order polynomial mean square fit coefficients and the zero-Hz 
Fourier component with transient related effects, the higher- 
order Fourier series terms have not previously been noted as 
particularly prominent in the raster signal display. In 
view of the series of nonlinearities present in obtaining the 
interpolated minimal phase raster points, it is not impossible 
that these higher order terms owe their separating capability 
to the data modulation of the signal source. In this respect 
the "sideband" structure of the single feature separation 
data of Figure 27 is most intriguing. 

Before completely discounting the usefulness of bias 
measures and transient phenomena, a more careful represen- 
tation and analysis of the data in the transient region should 
be attempted. The representation could perhaps involve the 
calculation of the minimal phase by use of a moving-average 
technique which v;orks backv/ard from data in the steady-state 
region to the starting time of the signal . At the same time 
use of a different decay constant and a smoother integration 
to obtain the Laguerre polynomial coefficients might also 
prove of value. 

If additional accurate ground truth data is available, 
the modeling of the probability distributions by something 
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other than a Gaussian form could be useful. The average 
squared magnitude of the Fourier components might be more 
accurately represented by a Rayleigh distribution. Although 
some of the measured parameters for one or more classes may 
not have a unimodal distribution, the number of accurately 
labeled samples currently available provides little encour- 
agement for the use of Parzen window probability estimators. 

The foregoing comments suggest that four or five parameters 
characterizing the raster scan data may provide signal iden- 
tifying information useful to about the same degree of accura- 
cy as visual interpretation of the raster scan display. Since 
the optimal parameters measured in this investigation do not 
appear to correspond to the visual cues used by system opera- 
tors, one has reason to believe that other parameters 
derived in subsequent analysis might improve upon this 
degree of accuracy. 

The measurement technique investigated exists outside the 
mainstream of the Parameter Encoder's indentif ication process. 
As such it constitutes a relatively slow procedure (one minute 
per signal) , requiring the use of disc data files which con- 
tain much of the same information available directly from 
the master file structure. Even should the programs used 
be modified to work solely with the master file, the data 
obtained might be of greatest use in a sequential identifica- 
tion scheme. It could be used much as the present raster 
scan display is; to resolve ambiguities in the choice between 
a relatively small number of classes. 
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APPENDIX A 



METHODS OF PROBABILITY ESTIMATION 



1 . Parametric Techniques 

Perhaps the most common technique which has been used 
to infer the a priori probability structure presupposes that 
the distribution of features in a given class is multivariate 
normal. As a first step the mean and variance for each feature 
of a class are estimated. For the class having samples 



variance vector represents a simplification from the more 
general multivariate form: 



of the d X d covariance matrix, which may be estimated by: 




A-1 




(i) ^2 



A" 2 




The 




The components of the vector are the diagonal elements 




A-4 
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Under the assumption of independent features, the 
distribution is characterized by the product of d Gaussian 



distributions . 



d 





A- 5 



If required, the technique may be extended to include the 
estimation of a covariance matrix for each class by use of 
Equations A-3 and A-4 . This procedure has enjoyed great 
popularity, particularly in the modeling of features with 
unimodal class conditional probability distributions. The 
ease with which the probability measure of a Gaussian distri- 
bution may be interpreted as a v;eighted distance measure 
has been an important factor in its retention. 

2 . Non-Parametr ic Techniques 



priori probability density function involves its direct 
estimation through the use of window functions or potential 
functions. In this technique, the estimate at some point x 
in the d-dimensional feature space is obtained by the super- 
position of terms which are generalized distance functions 
of that point and the training set points of the class 



Another technique technique for obtaining the a- 





A-6 
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The technique arises from a viewpoint proposed by 
Parzen. In one-diraensional cases, one can estimate the 
probability density at some point x by counting the number 
of samples N^^Cx) in the interval [x+h,x+h] and dividing 
by the product of the total number of samples M and the 
interval width 2h 



^ M2h 

If one defines a "window" function 



A- 7 



K(y) = 



7 - |y| 1 1 

0 |y| > 1 



A-8 



one can write 



Pm(x) = 



Mh 



M 

Z 

j=l 



K(- 



x-y. 



A- 9 



This corresponds to the potential function 



Y(x,y) = ^ K A-10 

Additionally, the argument presented may be extended to 
multiple dim.ensions by considering 
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K(u) 



1 , l^j I 1 § 3 = 1 , . . . ,m 

0 ^ otherwise 

1 2i “ y 

<2i/y) = ^ ) A(ll-12) 

'^m m 

where is a multidimensional volume of the hypercube of 

side h . 
m 

The window function K need not be the exact form of 

Equation A-11. It can be shown that if the estimate Pj^(x) 

is a function both of the total numb>er of samples M and 

V = V (M) , the conditions 
m m 

0 < K(y) < ” 



/K(y) dy = 1 
lim V = 0 

fk^-vco ^ 

lim MV = 

ensure that as the number of total samples M->-“ 

lim ^w(x) = P(x) 

lim (Pw{x) - p(x))^ = 0 

M-VCO 



A(13-16) 



A- 17 



A-18 
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m 

ij:-! 



Thus, we can consider Pj^(x) to be a "blurred" or "noisy" 
version of p ( x) as seen through an averaging v/indow. 

Intuitively, one can describe the following useful 
properties of the potential function: y(x,y) 

1. y(X/Y) should be maximum for x = y . 

2. y(2i»y) should be approximately zero for x distant 
from y . 

3. Y(x»y) should be continuous and decrease monotonically 
with the distance between y and x. 

4. If Y(x 3 ^»y) = where y is a sample point, the 

patterns represented by x^^ and X 2 should have about 
the same similarity to y. 

Such properties can be obtained with a potential function 
of the form: 

1 1 ^ ^^i " ^i^^ 

Y(2i'y) = ••■■■ (372 1 2 ^ 

2u a, ... a , 1=1 a • 

id 1 

where are arbitrarily chosen values. Although the form 
here is Gaussian, and the function of a distance measure, 
it is also possible to construct a potential function of 
the form: 

^ 2 

Y(x,y) = Z X/ 4). (x) (f). (y) A-20 

- i=l 1 1 ^ 

where is a complete set of multivariate orthonormal 
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functions and are constants. Such techniques have been 
described more completely in the work of Aizermann and 
Bravermann and others [1, 7, 19]. 

Given sufficient samples, the Parzen window or 
potential function approach essentially assures satisfactory 
convergence to an arbitrarily complex distribution. Unfor- 
tunately, this sufficient number of samples may be far greater 
than the number required if the form of the distributions 
was known. Since every sample point is used in the construc- 
tion of the density function, the above approach affords 
little economy in the way of data reduction and leads to 
a demand for computation time and storage space exponentially 
increasing with the number of features. 
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APPENDIX B 



CLASSIFICATION ALGORITHM 
1 . Distance Based Classification Schemes 

An assumption commonly used to implement a decision 
rule is that of multivariate normal class dependent density 
functions of the form; 



p(x|C.) = p(C.)— 37^ expt- itx-y^^h'^sT^[x-y^^h 

“ ^ ^ 2tt'^/^|S. ' 2 - - -1 - - 



1 



= p(C.) ^ exp t- i[x-il^'’h'^S^^(x-y^^h B-1 

i 



where 



X is the d-component feature vector 
y^ is the d-component class mean vector 
is a d X d class covariance matrix 
Is. I is the determinant of S.. . 

For such a distribution one may obtain the locus of 

points of constant density as a hyperellipsoid for which 

(i) T -1 (i^ 

the quadratic surface (x“]i ) (2£~iL ‘ constant. 

The quantity 



, (i).T^-l. (i). 

= (x-y. ) (2S"iL ) 



B-2 
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is sometimes called the Mahalanobis distance from x to . 

Since the probability densities are always non- 
negative, the maximum likelihood criteria 

p(C^)p(C^|x) ^ p (C^) p (Cj 1 x) for all j B-3 

implies 

log p(C^)p(C^Ix) > log p(C^)p(C^|x) B-4 

Substituting the multivariate density from Equation 
1, the above expres"sion yields 

log p(C^) + log - |[x-/]^s"^[x-/] 

- jji {log p(C^) - log V. - } 

B-5 

The above equation defines the decision boundaries 
of each class. If a common covariance matrix S^ = S^ = S 
is assumed for each class in the above equation, the log 
terms may be ignored and the form of the decision boundary 
becomes 



1. (i)>T„-l, (i), 

■ 2 -(x-y ) S (x-y ) 



P(C.) 



< min{ [x-y ^ ] *^S ^[x-y^^^] + log i— } B-6 






P(Cj) 
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Similarly, if equiprobable , a-priori class proba- 
bilities are assumed, p(C^) = p(Cj) then the decision rule 
becomes 

11) _ min{ [x-K^^hs"^Ix-y^ B-7 

This is equivalent to simply assigning the feature vector 
X to that class whose pattern, , is at the minimum 

Mahalanobis distance from x. The concept of Mahalanobis 
distance is particularly applicable to well-separated 
unimodal feature distributions, but presents a considerable 
problem in those cases where the distributions exhibit 
more than one local maximum. 

2 . K - Nearest Neighbor Algorithms 

The concept of distance may be applied in a slightly 
different form to probability functions not well characterized 
by the multivariate normal form. In this application, dis- 
tances in feature space are quite frequently normalized by 
some form of global variance or maximal spread. The nor- 
malized Euclidean distance from the sample point x to each 
point of the training set (x^ 

1 m 

calculated and ranked in order of increasing distance. In 
the simplest algorithm, class membership of the sample point 
is determined by the class membership of the greatest number 
of the k-nearest neighbors in the training set. For example, 
suppose that of the ten points of the training set nearest 
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to the test sample, three were members in class 1, four of 
class 2 and one each were members of classes 4, 5, 6, the 
simplest algorithm would choose class 2 as the probable 
classification. Alternative algorithms of this type weight 
the contribution of the k-nearest training points according 
to their distance ranking. The nearest point might be 
assigned a weighted class vote of 10, the next nearest, a 
weight of 9 and so forth. Class membership of the sample 
point would then be determined by summing the weighted class 
votes of the ten nearest neighbors and choosing that class 
which received the maximum number of votes. 

Such procedures are reminiscent of the Parzen window 
estimators described in Appendix A. 2, and in fact can be 
sho\nx to be equivalent if, instead of rank order, distance 
weighted contributions to the class vote are allowed. Numer- 
ous techniques [3, 4, 10] exist to increase the computational 
efficiency of these K-nearest neighbor algorithms by reducing 
the nvimber of Euclidean distance terms v?hich must be 
calculated for each sample point. It must be mentioned 
that these techniques are of greatest use when the individual 
classes are equiprobable or at the very least, when every 
class of the training set has more than k elements. 
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APPENDIX C 



FEATURE SELECTION AND RANKING 



1 . Linear Combinations of Features 

The objectives of low dimensionality and retention 
of sufficient information may often be obtained through the 
use of custom orthnormal transformations of the measurement 
space [13, 15, 19, 20], In such an interpretation the sample 
data is presumed to be adequately represented by the linear 
combination of some finite set of functions 



z 



n 

E ai _^i 

i=l ^ 



C-1 



where 



ii “ / 

For a given set of functions the weights are adjusted so 
that the average mean square error 



n 

E 

k=l 




n 

E 

i=l 



a . d) . , ) 



C-2 



is minimized. 

The values a^ are thus determined by taking partial 
derivatives with respect to a^ and setting the result equal 
to zero, obtaining n equations in n unknowns of the form 
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C-3 



n m 

E a. I 
i=l ^ k-1 



^ik *^jk 



m 



E z. 
k=l ^ 



4 > 



ik 



Now if the functions (ii . . . <b are orthonormal; that is if 

-*-i -^n 



ra 



Z (j) 
k=l 



ik 




0 

1 



i i- j 

i = j 



C-4 



Equation C-3 reduces to 



a 



i 





C-5 



Then values of a^ may then be used as the pattern space x 
representation of the vector z under the linear transformation 



X = T £ C-6 

T = [ '^2 ' * ' ’ '~n^ 

Expansion of a function by a Fourier series represents one 
form of orthonormal transformation, but one v;hich has 
distinct advantages in computational efficiency. 

It can be shown that Karhunen-Louve expansion of 
N m-dimensional samples z. = ) leads to a m x m 

eigenvalue problem. 

s Ik = 
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where 




C-8 




The above solution does in fact minimize mean square error 
averaged over the entire sample set. By choosing eigen- 
vectors corresponding to the n largest eigenvalues, one can 
achieve dimensionality reduction from m to n. This approach, 
also known as principal factor analysis, requires the compu- 
tation of the eigenvectors and eigenvalues of a large (m x m) 
matrix. Note that if the number of signals N is less than 
m, there will be at most (N-1) non-zero eigenvalues, and the 
minimum mean square error is zero. 



linearly independent patterns typical of the different classes 
or types of samples; then apply the Gram-Schmidt orthonormali- 
zation process to them directly. In this process one begins 
with the "typical" samples, perhaps derived from class means. 



A simpler procedure is to choose a small number of 



y 



... 

1 m 




y 




C-9 
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and constructs 



the orthonormal functions in the following 

manner : 

_1 
( 2 ) 

=12 = 'll'l ) 

*2 ' - =12 tl 

=jk ' 

_|Pk = y - c^k^i " 

_1 

ik “ <±k'ik> ^ *k 



*^k-l k -^k-1 

C-10 



ll = 
ii = 



where 



m 



(g.,g ) =— E g..,' 



'ki 



C-11 



Clearly, all of the samples used to generate the orthonormal 
vectors may be represented exactly by an n-term expansion; 
thus the representation error of the 'typical' sample is 
minimized. If these samples are truly typical, one can 
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expect other samples to be represented closely by the 
orthonormal basis functions. 

The dimensionality of the resultant pattern space 
is the same as the number of linearily independent typical 
samples . 

2 . Criteria for Feature Selection and Ranking 

Rather than achieving a pattern space through use 
of a transformation which minimizes mean square error over 
the training set, it is often desirable to optimize some 
measure of separation quality. Meisel [14] catalogs more 
than ten different measures of separation quality which are 
directly calculable from general a priori conditional density 
functions. For unimodal distributions, however, it is 
convenient to describe such a measure in terms of average 
distances in the feature space. 

Two distance measures need to be considered. One, 
the inter-class distance, might be characterized by either 
the distance beteen the means of the classes or the average 
squared distances between points of different classes. 
Another, the intra-class distance, is the measure of the 
internal scatter of samples and may be characterized by an 
average squared distance between the points of each class 
and summed over the classes. In the two-class problem, a 
common definition of the interclass distance, S, and intra- 
class distance, R^, is 
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s = 


, \ ‘'z' ’'s' xf 1 

^1 ^2 i=l j=l 


R. = 
1 


K. K. , 

1 1-1 3=1 k=3 ■’ 


2 

where d (x,y) 


is some distance measure of the vector 



variables x,y . . 

Conventional distance measures might use one of 



the following 


forms : 


Euclidian 


2 2 
d, (x,y_) = E (x. - y.) 

i=l 1 1 


city block 


n 

^2 = 2 Ix. - y . 1 

i=l 


Mahalanobis 


2 

d 3 (x»y) = (x - y) ~ (2£ " 


localized 

distance 


d^*(x,y) = 1 - exp[-i d^^(x,y)] C (14-17) 



In generalizing this distance criteria to a m-class problem, 
it is convenient to resort to the concept of scatter matrices 



where 




S. = 

~i 


Kf 

= E (x, - y^) (x, - y C-18 

k=l ^ ^ 


§w = 


m 

= E S. C-19 

i=l''^ 
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!: A 

* 




represent the individual class scatter matrices and the 
collective within-class scatter matrices respectively. If 
one defines a total scatter matrix in terms of a global 
mean M 



M 



I 

i 



K 



m 

I 

i=l 



K. y 



(i) 



C-20 



S = E (x - M) (x - 
all samples 



C-21 



we find that 



S„ = 



m 

E 

i=l 



E (x - 
xeC^ 



) (x - 



(i) ,T 

y ) 



m 

E 

i=l 



E(y 

xeC^ 



- M) 



-w 



+ 



:b 



C-22 



provides a natural definition for both inter-set and intra- 
set scatter matrices: 



m 

S = E S . 

.W 1 



C-23 



(i) i T 

= E K, - M) - M) ^ 



i=i 1 



C-24 



These three matrices not only provide the basis for multiple 
discriminant analysis [7, 14, 18], but their trace or 
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determinant can also serve as a measure of separability for 
use in an algorithm involving either linear combinations 



of features or discrete feature selection. For example, 
as criteria for separability, one might choose to maximize 
the quality function 



Q 



1 



tr S 
~w 



C-25 



or to minimize one of the quality functions 



Q 



2 






-1 



S 

~w 



C-26 



Q 



3 




C-27 



All of the above criteria are particularly interesting in 
that they can be evaluated as a result of multiple discrim- 
inant analysis using linear combinations of features. 

For a c-class problem we can accomplish the pro- 
jection of the d-dimensional feature space onto some alterna- 
tive space of not more than c-1 dimensions through the use 
of a non-unique d x (c-1) matrix W. 



y = 




C-28 



It is straightforward problem to show that the resultant 

scatter matrices S , S in this projection space are 

^ W -V 13 
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C-29 



S = W S W 






If one then chooses W to maximize the ratio of between-to 
within-class scatter matrix determinants, 



«4 ' 



I5bI 

lii 



!b 

|w'^ w| 
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there is a solution of the form 



5 b Hi 



= X. 



S„ w. 

,vW — 1 



C-31 



If is non-singular one can directly solve an eigenvalue 
problem related in form to the Karhunen-Louve problem 
Equation C-7. 




;b 



w 

— 1 



X. w. 

X —X 



C-32 



or one may then solve the characteristic polynomial 



lS„ - X. S I = 0 C-33 

' ~B 1 

for the largest eigenvalues and find the corresponding 
eigenvectors from 
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0 



C-34 




X. 

1 



s 






w 
— 1 



If m-1 of the vectors - M) are linearly 

independent, there will be not more than m-1 non-zero eigen- 
values for this problem. If the within-class scatter 
matrix is isotropic, which is the case if all features are 
independent and suitably normalized, the resultant eigen- 
vectors span the space defined by the vector set - M) 

and may be established by Gram-Schmidt orthogonalization . 

The eigenvalues of this problem serve to characterize the 
trace and determinant criteria. For example, the use of 
all features leads to the criteria values 
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As was previously mentioned, the criteria described 
here are best suited to unimodal distributions and increase 
the theoretically achievable error rate. In those cases 
where it can be established that a unimodal distribution is 
inappropriate, the use of separability measures more closely 
related to the estimated conditional probability density is 
indicated [16] . For example one might use a quality measure: 
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i=l ^ ^ 



where N is the total number of training set samples and 
F(x) represents some selection of features from the initial 
feature space and 

m 

P(x) = I p(C.) p(xlc.) C-39 

i=l 

Qg is in fact a measure of overlap and has an optimal value 
of zero. 

3 . Feature Search Algorithms 

A cursory inspection of the problem of finding those 
n of m features which optimize one of the criteria given in 
the previous section shows that 
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evaluations must be performed if one is to conduct an 
exhaustive search of all possible feature combinations. The 
exhaustive search for the best ten or twenty features for 
example, requires consideration of more than 184,000 possibil- 
ities. Typically, the practical feature search employs 
siiboptimal procedures that may be justified if simplifying 
assumptions regarding the nature of the feature distributions 
are allowed. 
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The simplest suboptimal search procedure is to assume 
the independence of features and evaluate the desired 
criteria function for all single features. These m features 
are then ranked in order of decreasing optimality and the 
first n are chosen. The number of criteria evaluations is 
the same as the dimension of initial feature space. 

A second search method used by Mucciardi and Gose 
[17] involves a technique known as search without replacement. 
This method requires N^, criteria evaluations where 



M _ n (m - n) 
E 2 
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Fu [9] has proposed a technique of sequential feature 
selection which can be modified to include the cost of 
feature measurement and Chang [3] has developed alternative 
dynamic programming approaches to feature selection, one of 
which requires only a slightly greater number of evaluations 
Np, than the search without replacement algorithm. 



Ng = n (m - — 2 ^ ') + n - 2 (m - 1) 



C-42) 



The single feature ranking technique has the 
disadvantage of ignoring the effects of feature correlations. 
The search without replacement technique provides a method 
for treating these correlations, but it also presumes that 
the features obtained through n stages of conditional single 
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feature evaluation procedure are the n best features, which 
need not be the case. In general, as one progresses to 
successively more complicated search procedures, more of 
the effects of feature dependence may be taken into account. 
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PARAMETER HISTOGRAMS 

o 






C\J 

CM 

CM 



rA 



CO tN-MDM:) LP\ 
OOD 

CT\ CTn O O O- D'- L^^ 
Cn £>- O- O- LA rA 

0^ C^ CTn LA nA ^A NA 

A-hA 
NANAiA NA 
OOOOOO rAiArANA 
QO 
CO 

JA 



00 






CM 

O 

C5^ 



a 

O 

•H 

-P 

0 

rH 

o 

W 

Q) 



e 

fd 

p 

cn 

o 

-p 

w 



ro 

ro 

s 

D 

O 

H 

p4 



1 

LO 



130 



LTv 



V- O 



GO 

OJ 



lA 



CO 
LA 
■ (A- 






O.I>- 

!>- 
rA 
OnA- 
0^ A- 
O^A-LA 
A~ 

COCO 
C^^ 
rArA 
C^rA 
rA 
N^rA 
rHrANMA 
r^rvj 



LO 



r-lC^ rA 
VX) C\J 



OJ 

VD 



VO 



VO 



CO 



VO 



CO 

CO 



o 

LT> 



LO 

I 



o 

I 



131 



FIGURE 34. Histogram: Zero Order Coefficient of Polynomial Fit (903 
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FIGURE 35. Histogram: 0 Hz Frequency Component (918) 



o 

CO 



LO 

CM 



O 

CM 



CU 




0-[>- 



LA 

C\J 

O-A- 

VO 



vx>rvj 

A- 

LA 



o 



CTn 



CO « — < — 
A-LAhA 



CO 



CO ^ rA 
CO rArArA 
CO A-LArA 
C^ 

CO CAA-M3 
O^C5^^A^A^A 
CTnvO rA 
CArA 
00*“ 
O^K^^A 
VO hA 
A-hA 



1 

LO 



o 



-L- 

LO 



CA 



o^ 

—7 



LO 

CM 



133 



FIGURE 36. Histogram: 18 Hz Frequency Component (926) 
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FIGURE 37. Histogram: 20 Hz Frequency Component (927) 
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FIGURE 38. Histogram; 28 Hz Frequency Component (931) 



COMPUTER PROGRAMS 



THE PROGRAMS DISCUSSED IM THIS SECTION ALL USE 
LIBRARY AND EXECUTIVE ROUTINES WHICH ARE PART OF THE PARA- 
METER ENCODER SOFTWARE PACKAGE. AS AN AID TO INTERPRETATION 
A BRIEF DESCRIPTION OF SOME OF THESE ROUTINES FOLLOWS 



ROUTINE FUNCTION 

ALPHA CONTROLS TYPE OF OUTPUT TO CRT 

CHOUTdl OUTPUTS SINGLE CHARACTER TO CRT BUFFER 

CURSRI IC, IX, lY) 

DISPLAY CURSOR ON CRT 

HALT UNTIL KEYBOARD DEPRESSED 

IC= KEYBOARD CHARACTER 

IX= X CURSOR POSITION 

IY= Y CURSOR POSITON 

OSK2( 14, NBUF,NWD, JR) 

14/15 = READ/WRITE MASTER FILE 
NBUF= BUFFER NAME 
NWD= NUMBER OF WORDS 
JR= RELATIVE SECTOR OF MASTER FILE 

ERASE CLEARS CRT DISPLAY 

EXEC(10,NM) CAlL PROGRAM STORED IN DISC FILE NM 

EXEC( 14, IDCON,NBUF,NWD,NAME,NDX) 

EX EC ( 1 5 , I DCON , NBUF , NWD , N AME , NOX ) 

14/15= READ/WRITE DISC FILE 
IDCON= DISC CONTROL WORD 
N8UF= START LOCATION OF BUFFER 
NWD= NUMBER OF WORDS IN BUFFER 
NAME= NAME OF FILE TO READ/WRITE 
NDX= STARTING SECTOR OF FILE 

JFINDI 1ST, I END, IV AL, I FID 

CONDUCTS BINARY SEARCH FROM 
IFIL(IST) TO IFIL(IEND) 

RETURNS INDEX OF LAST VALUE 
LESS THAN OR EQUAL TO IVAL 

JPLOT (K, IX, lY) 

PLOT ON CRT 

K = CONTROL OF POINT OR LINE PLOT 
IX= X POSITION 
IY= Y POSITION 

POSAC(K) POSITIONS ALPHANUMERIC CURSOR 

K= LINE NUMBER 

SL2RS( JMFO,NS, JS, JR) 

FIND RELATIVE SECTOR FROM MASTER 

FILE DIRECTORY LOCATION 

JMFO= MASTER FILE INFQRMATICN BUFFER 

NS= NUMBER OF SIGNALS TO FIND 

JS= DIRECTORY LOCATION 

JR= RELATIVE SECTOR IN MASTER FILE 
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interface program 



FTN 

PROGRAM FACE 
DIMENSION IFIL(3) 

DIMENSION IFL(3) 

DIMENSION NM(3) ,NAME(3) 

DIMENSION NAM(3) 

COMMON M,MDX, ISIG( 126) 

COMMON RATE, N, IX( 256) 

COMMON IC(30) 

COMMON IDR(768) 

EQUIVALENCE ( T I MB , IC ( 7 ) ) , { TT , I C I 12 ) ) 

EQUIVALENCE (TOT, IC(14)),(RST,IC(16)) 

EQUIVALENCE (DIS,IC(20)) 

THIS PROGRAM REFORMATS THE DISC FILE PARF TO 
A FORM ACCEPTABLE BY THE ESL PROGRAM RASF 
IT SETS UP A FILE PMAP TO PASS INFORMATION TO THE 
RASTER DISPLAY PROGRAM 

VARIABLE INFORMATION 

IFL NAME OF FILE CONTAINING FORMATTED INFO 
NAM NAME OF FILE TO BE REFORMATTED 
NM NAME OF NEXT PROGRAM 
NAME NAME OF CURRENT PROGRAM 

RATE, N, IX BUFFER FOR SIGNAL ENTRY FROM PARF 
RASTER PERIOD, NUMBER OF ENTRIES, 

DELTA TRANSITION TIMES 
IC(l) NUMBER OF WORDS OF FILE INFO 
IC(2-4) NAME OF FILE ^^,'HERE TIME DATA STORED 
IC(5) # SECTORS IN FILE PMAP 

IC(6) # SIGNALS IN FILE PMAP 

IC(7-8) QUANTIZATION INTERVAL 
IC(9-ll) FILE TYPE CODING 
IC( 12-13) SIGNAL DURATION 

IC(14-15) NUMBER OF TRANSITION TIMES IN FILE 

IC(I6-17) NOMINAL RASTER PERIOD 
IDR DIRECTORY BUFFER FOR PARF 

DATA IFL( 1 ) , IFL( 2) , IFL(3) 

X /50115B,40520B,20040B/ 

PM AP 

DATA NAM( 1) ,NAM( 2) ,NAM(3) 

X /50101B,51 106B, 20040B/ 

PA RF 

DATA NM( 1 ) ,NM(2 ) ,NM(3 ) 

X /51 101B,51506B,20040B/ 

RA SF 

DATA NAME! 1) ,NAME( 2) , NAME( 3 ) 

X /43101B, 41505B, 200403/ 

FA CE 

CALL EXEC ( 14,102B,IDR,768,NAM,0) 

1 CONTINUE 

READ(5,200) ISG 

200 FORMAT ( 14) 

IF(ISG) 1,999,201 

201 CONTINUE 
DO 10 1=1,768 
IF ( IDR ( I)- ISG) 10,5,10 
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5 NDX=I*2+4 
GO TO 20 
10 CONTINUE 

WRITE (6,101) ISG 

101 FORMAT ("SIGNAL NO. ",I5,"N0T PRESENT") 
GO TO I 



20 CONTINUE 

CALL EXEC ( 14,102B,RATE,256,NAM,N0X) 
CALL EXEC (15,102B,IX,256,IFL,0) 

CALCULATE SIGNAL DURATION 

TT=0. 

DO 45 1=1, N 

IF ( IX( I ) ) 45,45,44 

44 CONTINUE 

TT = TT<-FLOAT( IX( I ) ) 

45 CONTINUE 

FORMAT FILE INFORMATION 

IC( 1)=21 
DO 40 1=1,3 
K=I+1 

IC(K)=IFL( I ) 

40 CONTINUE 
IC( 5)=2 
IC(6)^1 
IC(9)=52117B 
IC( 10)=0 
IC( 11)=1002B 
IC( 18)=0 
IC( 19)= ISG 
RST=-RATE 
TIMB=.0001 
RST=RST*TIMS 
TOT=FLOAT(N) 

DIS=-1. 

SET UP FILE PMAP USING ESL ROUTINE 

CALL USENT ( IC,NAME) 

CALL EXEC (10, NM) 

999 STOP 
END 
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MEASUREMENT PROGRAM 



FTN 

PROGRAM FOURD 
DIMENSION FTX(256) 

DIMENSION NAMRI3) 

DIMENSION NPAR(3 ) 

COMMON 10(2048) 

COMMON BP (236) 

COMMON SIG(128) 

COMMON I0M(3) .SNUMtIOMYdl) 

COMMON NAME(3) 

COMMON TT,T1,T2 
COMMON RATE, PER, TBtMOOE 
COMMON IDMM(9) 

COMMON JX,JY 
COMMON MODI 
COMMON ISG 

EQUIVALENCE ( FTX ( I ) , 1 0 ( 1025 ) ) 

EQUIVALENCE {NDX,BP(256)) 

THIS PROGRAM IS A DRIVING ROUTINE WHICH 
CALLS THE MEASUREMENT AND INTERPOLATION 
SUBROUTINE ,FAST FOURIER TRANS FORM , AND 
THE FOURIER SPECTRUM DISPLAY 

THE FOLLOWING VARIABLES ARE USED: 

FTX COMPLEX INPUT TO FAST FOURIER TRANSFORM 
SIG MEASURED SIGNAL PARAMETERS 

DATA NAMR{ 1) ,NAMR(2) ,NAMR(3) 

X /43101B ,41505B ,200403/ 

FA CE 

DATA NPAR( 1) ,NPAR(2) ,NPAR(3) 

X /50101B, 510408,200408/ 

PA R 

ZERO PARAMETER MEASUREMENTS 

DO 5 1 = 1, 128 
5 SIG(I) = 0. 

CALL RESDU 

PAUSE FOR RASTER DISPLAY 
PRESS Q TO QUIT 

CALL CURSR( iC, IX, lY) 

IF(IC-121B) 10,20,10 
10 CONTINUE 

CALCULATE 128 COMPLEX POINT FAST FOURIER TRANSFORM 

CALL F0UR2(FTX, 128,1 ,-l ) 

CONTINUE 
CALL FORD 

PAUSE FOR FOURIER DISPLAY 

CALL CURSRI IC, IX, lY) 

IPdC-dlB) 30,20,30 

WRITE PARAMETERS TO DISC FILE PAR 

30 CALL EXEC( 15, 102B ,SIG,256,NPAR ,NDX) 
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CALL PROGRAM FACE FROM DISC 

20 CONTINUE 

CALL EXEC (lOjNAMRJ 
EN0$ 



SUBROUTINE RESDU 
DIMENSION IT(256) 

DIMENSION TX(2J 

DIMENSION FTX(256) ,XF(2,4) ,B(5) 

DIMENSION X( 2,128) 

COMMON 10(2048) 

COMMON BP(256) 

COMMON SIG(128) 

COMMON IDM(3),SNUM,I0MY(li) 

COMMON NAME(3) 

COMMON TT,T1,T2 

COMMON RATE, PER, TB, MODE 

COMMON LCT,SCX,SCR,STX,STR 

COMMON JX,JY 

COMMON MODI 

COMMON ISG 

EQUIVALENCE ( X ( 1 , 1 ) , I D ( 1 5 37 ) ) , ( F TX ( 1 ) , I D ( 10 25 ) ) 
ECUIVALENCE( IT, ID( 513) ) 

EQUIVALENCE(TX,T1) 

POLY(F) =(SIG(5) 4-F=!MSIG(6) fF«(SIG(7) 
1+F>!‘(SIG(8)+F«SIG(9) ) ) ) )/PER 

THIS SUBROUTINE INTERPOLATES BETWEEN PTS OF THE 
SCAN, PROVIDES MEAN SQUARE FITS FOR EACH BIAS 
AND OBTAINS A RASTER VARIANCE MEASURE 

THE FOLLOWING VARIABLES APPEAR IN COMMON 
ID(l-256) delta TRANSITION TIMES 
IT TRANSITION TIMES 

X(1,K) RASTER PHASE POINTS 
X(2,K) FLOATING PT TRANSITION TIME 
SIG SIGNAL PARAMETERS MEASURED 
SNUM TOTAL NUMBER OF TRANSITIONS 
TT SIGNAL DURATION 

RATE , PER, TB, MODE, LCT, sex, SCR ,STX,STR,JX,JY 
CONTROL VARIABLES FOR DISPLAY 
ISG SIGNAL IDENTIFIER 

XF(2,4) 4 POINTS USED FOR CUBIC INTERPOLATION 

B COEFFICIENTS OF INTERPOLATED FIT 

INITIAL IZATION 

NTOP=IFIX( SNUM) 

LT=0 

RX=JX 

RY=JY 

LX=900 

LY=500 

SCAX=FLOAT(LX) /TT 
RAT2=PER/2 . 

RAT4=RAT2/2. 

CONVERT DELTA TIME TO TIME 

00 200 I=1,NT0P 
IF(IO(D) 200,200,201 
201 I0EL=I0(I) 

LT=LT+-IOtL 
200 IT(I)=LT 

FIND STARTING INDEX FOR MINIMAL PHASE DIFFERENCE 
JST=JF I ND( 1 ,NTOP ,400, I T) 
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S IG(2) =FLOAT( I T ( NTOP -2 ) -400 ) 

DELT= SIG(2)/128. 

SIG(2)=1./(SIG(2)*TB) 

N0X = 1 

NSWCH=1 

OHSQ=0. 

JT0P=0 

KST=1 

999 CONTINUE 
PHM=0. 

PHSQ=0. 

C 

C INTIALIZE RASTER PHASE POINT 

C 

J = JST 
210 J=J-2 

IF(J-l) 212,212,209 

212 PHH=0. 

TST=0. 

GO TO 214 

209 IF(ID(JI) 210,210,213 

213 TST=IT( J) 

PHH=AMOD( TST, PER) 

IF(PHH-RAT2) 214,215 

215 PHH=PHH-PER 

214 X(1,1)=TST 
X(2,1J=PHH 
DELL=0. 

PHM=PHM+PHH 

PHSQ=PHSQ+PHH*PHH 

N=TST/PER 

J = 1 

PHL=PHH 

C 

C CALCULATE MINIMAL PHASE DIFFERENCE POINTS 

C 

DO 300 I=JST,NT0P,2 
IFdDdl) 300,300,302 
302 TST=IT( I) 

41 N=N-H 

PHH=TST-FLOAT(N)*PER 

DEL=PHH-PHL 

IF(ABS{ DED-RAT2) 42,42,41 

42 N=N-1 
C 

C UPDATE MEANS AND VARIANCE ESTIMATES 

C 

PHM=PHM<-PHH 
PHSG=PHSQ+ PHH*PHH 
OHSQ=OHSQ^■DEL=i-'DEL 
IF(A6S(D£L)-RAT4) 305,707 
707 DEL=DELL*.5 
PHH=PHL+DEL 
305 J=J+1 

X( 1, J)=TST 
X(2, J)=PHH 
DELL=DEL 
PHL=PHH 
300 CONTINUE 
C 

c 

c CALCULATE MEAN AND VARIANCE FOR BIAS CONDITION 

C 

JTOP = JTCP<- J 
JP = J 
XNUM=JP 
XNUM1= JP-1 

SIG(35) = SIG(35 )4-(PHSQ-PHM«PHM/XNUM)/XNUM1 
PHM1= PHM/XNUM 

c 

c 

c CALCULATE POLYNOMIAL MEAN SQUARE FIT 
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CALL CURV(JP) 

INTERPOLATE BETWEEN POINTS 
K = 1 

CALL CURF(X( ItK) ,3,B) 

XP=400. 

DC 601 1=1,128 
607 IF(XP-X(l,K+2n 602,603 

602 CONTINUE 
DX=XP-X( 1 ,K) 

FTXIKST ) = ( B( 1 ) +DX*(B( 2 )+DX«(B(3 ) +DX«B(4 ) ) ) ) /PER 
IX=XP=(‘SCAX + RX 

I Y=(FTX(KST)+.5I *FL0AT( LY) +RY 

DISPLAY INTERPOLATION 

CALL JPLOT (0, IX-4, IY-6) 

CALL ALPHA 
CALL CHOUT(52BJ 

INPUT FOR FAST FOURIER TRANSFORM 

FTX(KST)=FTX(KST )-P0LY(XP) 

KST=KST+2 
XP=XP+D£LT 
GO TO 601 

603 IF(K+4-JP) 604,602,602 

604 CONTINUE 
K = K^-l 

CALL CURF(X( I,K) ,3,B) 

GO TO 607 
601 CONTINUE 

CALL CH0UT(-1 ) 

PLOT POLYNOMIAL MEAN SQUARE FIT DATA 

XP=400. 

DO 20 1=1, 128 
1X=XP*SCAX+RX 

IY=(P0LY(XP) ^-.5) «FL0AT(LY) + RY 
CALL JPLOTIO, IX-4, IY-6J 
CALL ALPHA 
CALL CH0UT(53B) 

XP=XP+DELT 
20 CONTINUE 

CALL CHOUT(-l) 

CHECK TO SEE IF BOTH BIAS CONDITIONS ARE PLOTTED 

NSWCH=-NSWCH 
IF(NSWCH) 700,701 

700 CONTINUE 
XP=400. 

PHM2=PHM1 
JST=JST+1 
KST = 2 

S IG(4) = SIG(3) 

DO 702 1=1,5 

Kl=I+4 

K2=I+9 

SIGIK2 )=SIG(K1) 

702 CONTINUE 
GO TO 999 

701 CONTINUE 

CALCULATE MEAN BIAS AND INTRINSIC VARIANCE 
XNM=JT0P«2 

SIG(20)= ABS( PHM2-PHM 1 ) 

SIG(36)= OHSQ/XNM 



142 



oooooooooo oooo 



CALCULATE LAGUEKRE COEFFICIENTS 

P = 80. 

TEN0=800. 

CALL LAGER! ID, PER, P, TEND, SI G( 15 ) ) 

RETURN 

END 



SUBROUTINE CURV(NPTS) 

DIMENSION X(2,128) 

DIMENSION B(5) 

DIMENSION A(5,6) 

COMMON ID(2048) 

COMMON BPI256) 

COMMON SIG(128) 

COMMON IDMMY(16) 

COMMON NAME! 3) 

COMMON TT,T1,T2 
COMMON RATE , PER, TB , MODE 
COMMON LCT,SCX,SCR,STX.,STR 
COMMON JX JY 
CCMMON MODltlSG 

EQUIVALENCE (X( I , 1 ) , ID( 1537 )),(B(l),SIG(5)) 

THIS SUBROUTINE CALCULATES THE COEFFICIENTS 
OF A POLYNOMIAL MEAN SQUARE FIT OF ORDER LESS 
THAN FOUR, USING UP TO 128 SAMPLE POINTS 

X(2,128) SAMPLE POINTS 
SIGI5-10) POLYNOMIAL COEFFICIENTS 

SPECIFY ORDER OF POLYNOMIAL TO BE FIT 

NUM=2 

DO 20 J=l,6 
DO 20 1=1,5 
A(I,J)= 0. 

20 CONTINUE 

A( 1,2) =FLOAT(NPTS) 

DO 100 I=1,NPTS 
YT= X(2, I ) 

X1=X( 1,1) 

X2= Xl«Xl 
X3= X2*X1 
X4= X3«X1 
X5= X4*X1 
X6= X5*X1 
X7= X6*X1 
X8= X7*X1 
C 

A(l,l)= A( I,1)+YT 
A( 1,3)= A( 1,3)+X1 
A( 1,4) =A( 1,4)<-X2 
A( 1,5)=A( 1 ,5)+X3 
A(1,6)=A( 1,6)+X4 
A(2,1)=A( 2,1 )<-Xl«YT 
A(2,6)=A( 2,6H-X5 
A(3,1) = A(3,1 )+X2^YT 
A(3,6)=A(3,6)«-X6 
A(4,l)=A(4,l)+X3=!‘YT 
A(4,6)=A(4,6)+X7 
A(5,1)=A(5,1)<-X4«YT 
A( 5,6) =A( 5,6)+Xa 
100 CONTINUE 

DO 25 1=2,5 
DO 25 J=2,5 
A(I,J)=A( I-l,Ji-l) 
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25 CONT INUE 
LCT= NUM+2 
LMT= NUM+1 
DO 35 l=lfNUM 
DO 35 J = 1 , I 
L= J+1 
K = H-1-J 
M = K+1 

FAC= A(K,L)/A{M,L) 

CALL ROW( A,NUM,M,FACf K) 

35 CONTINUE 
M= LCT 

DO 45 1=1 »LMT 
00 36 K=1,LCT 
A( I ,K)=A< I ,KI/A( I ,M) 

36 CONTINUE 
IF(I-LNT) 37,45,45 

37 CONTINUE 

DO 38 J=I,NUM 
MK=J+1 
FAC=A( MK,M) 

CALL R0W( A,NUM, I ,FAC,HK) 

38 CONTINUE 
M= M-1 

45 CONTINUE 
J = LCT 

DO 46 1=1, LMT 
J = J-1 

6(JJ= A( I, 1) 

46 CONTINUE 
IF(LCT-5) 70,70,75 

70 DO 47 I=LCT,5 
B( I )= 0. 

47 CONTINUE 
75 RETURN 

END 



SUBROUTINE ROW ( A , NUM , I RM , FAC , I RS ) 
DIMENSION A(5,6) 

DIMENSION TEMP(6) 

K J=NUM^-2 

DO 50 1=1, KJ 

TEMPI I ) = FAC*A(I RM, I ) 

50 CONTINUE 

DO 100 1=1, KJ 

A( IRS, I )= A( IRS, I)-TEMP(I ) 

100 CONTINUE 
RETURN 
END 
ENDS 



SUBROUTINE LAGER ( I D, RATE, P, TEND SL) 

DIMENSION ID(256),SL(5) 

THIS SUBROUTINE CALCULATES THE LAGUERRE POLYNOMIAL 
COEFFICIENTS FOR THE FIRST FIVE POLYNOMIALS USING 
TRAPEZOIDAL INTEGRATION 

DO 5 1= 1,5 
5 SL( I ) =0. 

CALCULATE MINIMAL PHASE DIFFERENCE RASTER POINTS 
USING BOTH BIAS CONDITIONS 

PHL=0. 

TT= 0. 

N= 0 
P= l./P 
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RAT2 =RATE/2. 

DO 10 1=1,256 
IF(ID(I)) 10,10,11 

11 TT=TT + FL.OAT( ID( I ) ) 

IF(TT-TEND) 12,12,13 

12 DEL= FLOAK ID( I ) ) 

14 N=N+1 

PHH=(TT-FL0AT(N1*RATE) 

01F= PHH-PHL 

IF(ABS( DIF )-RAT2J 15,15,14 

15 N=N-1 

CALCULATE CONTRIBUTIONS TO LAGUERRE POLYNOMIALS 
K=1 

T1=TT-0£L 

30 K=K-1 

0 IF = PHL 
RT=P*T1 

REXP=.5*DEL-EXP( -RT) 

SL(1)= SL( l) + REXP=f=DIF 
SL(2)=SL(21 + (2.*RT-1.)=:'REXP={=DIF 
SL(3)=SL(3)<-( l.+RT*(-4.+RT*(2. ) I )*REXP’«=DIF 
SL(4) = SL(4)f(-i. +RT=M6. +RT* (-6. +R T’^ ( 4 . / 3 . ) ) ) ) *REXP*0IF 
SL(5) = SL(5) + ( 1 .+RT'-:M-8.+RT>:M 1 2 » +RT*{ -16 . / 3. +RT«2 . / 3 . ) ) 
1 ) )*REXP*DIF 

1 *REXP*DIF 
T1 = TT 
PHL=PHH 

IF(K) 31,30,30 

31 CONTINUE 
10 CONTINUE 

NORMALIZE 

13 PN0RM=SGRT(2./P) 

DO 20 1=1,5 
20 SL( I )=SL( I I^PNORM 
RETURN 
END 



SUBROUTINE FORD 
DIMENSION FX(256) 

DIMENSION NM(3) 

DIMENSION FTX(256) 

COMMON ID (2048) 

COMMON BP (256) 

COMMON SIG(128) 

COMMON IDM(3),SNUM,IDMY(11) 

COMMON NAME(3) 

COMMON TT,T1,T2 
COMMON RAT E,PER,TB, MODE 
COMMON IDMM(9) 

COMMON JX,JY 
CCMMCN MODI 
COMMON ISG 

EQUIVALENCE ( FX ( I ) , I D ( 1 537 ) ) 

EQUIVALENCE ( FTX (1 ) , I D (1025 ) ) 

EQUIVALENCE ( NDX , B P ( 256 ) ) 

THIS SUBROUTINE CALCULATES AND PLOTS FOURIER 
MAGNITUDES FOR THE 64 TERMS 

RX = JX 
RY = JY 
LX= 900 
LY= 500 
XMAX= 200. 

YMAX=4. 
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XMX4= XMAX/4. 
yMX2= YMAX/2t 
SIG(l)= RATE 
STX=0. 

STY=0. 

SCX= FLOAT(LX) /XMAX 
SCY= FLOAT(LY)/YMAX 
CALL ERASE 
CALL P0SAC(2) 

WRITE(6,100) ISG,SIG(37) ,SIG(2) 

100 FORMAT("FREQ MAG PLOT FOR SIG NO. '*,15/ 

1, "VARIANCE OF RE S I DUE = " t F 1 0 . 6 . 

2 " RESOLUTION= ",F10.6) 

CALL AXIS! JX.JY, sex, SCY.STX, STY, 

1 XMX4,YMX2,XMAX, YMAX) 

CALCULATE FIRST 64 MAGNITUDE COMPONENTS 

DO 40 1=1,63 

J=2*I+1 

K=257-2#I 

L = I + 1 

JI=J+1 

KI=K+1 

FX(L) = ( SQRT(FTX( J ) *>;^2 + FTX( K ) ={=*2<-FTX( JI )**2 
1 +FTX( KI )-^=2)J /2 . 

40 CONTINUE 

FX( 1 ) = ( SQRT(2.*( FTX( 1) >^=<<2 + FTX( 2)»«2) ) )/2. 
PLOT SPECTRUM 
DO 50 1=1,64 

I X=FLOAT( I-l )*S IG( 2J={=SCX + RX 

I Y=FX( I )*SCY+RY 

CALL JPL0T(1-I,IX,IY) 

50 CONTINUE 
CALL ALPHA 
CALL CH0UT(-1 ) 



DO 10 1=1,64 
K=I+39 

SIG<K)=FX( I) 
10 CONTINUE 
RETURN 
END 
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STATISTICS AND SCALING PROGRAM 



FTN 

PROGRAM CSTAT 
DIMENSION SPAR( 50) , JW( 50) 

01 MENS I CN NPAR(3 ) fNCLSO ) ,NL IB( 3) 

COMMON JSAVE(50) 

COMMON JMFU (128) ,JMFD(768)t MFE(334) 

COMMON IP0(256) , SIG( 128) 

COMMON ICD(128) , ISG(128) 

COMMON JBUF (4) ,CLMN(50 ) , CLVA(50 ) , JNPA( 50 ) , JPNUM(50)i 
1 JCD0(80) 

COMMON JGL0B<4) , GM ( 50 ) , GV( 5 0 ) , G TV( 50 ) , J BNO ( 50 > , 

1 JGDM(30) 

EQUIVALENCE ( SPA R ( 1 ) , MF E ( 5 ) ) , ( JWU) ,MFE( 105) ) 
EQUIVALENCE ( I D.NO , JBUF ( 1 ) ) , ( NEN T , J BUF { 3 ) ) 

This PROGRAM PLACES MEASURED SIGNAL PARAMETERS IN THE 
SIGNAL FILES OF THE MASTER FILE. CLASS STATISTICS 
ARE CALCULATED AND STORED IN THE MASTER FILE CLASS 
ENTRIES. 

THE FOLLOWING VARIABLES ARE USED: 

JSAVE STORAGE FOR OLD PARAMETERS 
JMFO INFORMATION ON MASTER FILE STRUCTURE 

JMFD MASTER FILE DIRECTORY 

MFE BUFFER FOR MASTER FILE ENTRY 
IPD PAR FILE DIRECTORY 
SIG BUFFER FOR PAR FILE ENTRY 

ICD CSIGF FILE DIRECTORY 

ISG BUFFER FOR SET OF SIGNAL NUMBERS OF SAME 
CLASS 

JBUF-JCDO BUFFER FOR ENTRY TO CLASS FILE 
JGL06-JGDM BUFFER FOR ENTRY TO GLOBAL FILE 



DATA NPARd > ,NPAR(2) f NPARI3 > 

1 /50101B, 51040B, 20040B/ 

PA R 

DATA NCLSI i> ,NCLS(2) ,NCLS(3) 

1 /41523B, 44507B, 43040B/ 

CS IG F 

DATA NLIBI 1) ,NLIB(2),NLIB(3) 

1 /46501B, 515248,430403/ 

MA ST F 

READ IN DIRECTORIES FOR FILE PA R,CS I GF , MASTF 

CALL EXEC( 14, 102B , IPD ,2 56, NPAR,0) 

CALL EXECI 14, 102B, ICD,128,NCLS,0) 

CALL EXECI 14, 102 B, JMFO, 123, NL IB, 0) 

NWD=JMF0(2 ) 

CALL DSK2( 14, JMFD,NWD, 1 ) 

CALL SL2RS(JMF0,1 , 1 ,LG) 

CALL DSK2I 14, JGL0B,384,LG) 

ZERO GLOBAL MEANS AND STND DEV 

DO 10 1=1,50 
GMI I )=0. 

GVI I )=0. 

GTVI I )=0. 

10 CONTINUE 
JGLOBI 3)=0 
CALL ERASE 
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CALL POSAC(O) 

WR I T E ( 6 , 1 1 3 ) 

113 FORMATC "ENTER INPUT DEVICE (5 FOR CR, 6 FOR TEK) ") 



READd,*) LU 
READILU,#) N 

CHANGE PARAMETER NAMES IN MASTER FILE INFO SECTOR 

DO 90 I 1 = 1 ,N 

Kl= 11+14 

JSAVEI I U = JMF0(K1 ) 

READILU;*) JMFO(Kl) 

90 CONTINUE 
JMFC( 14)=N 
1 CONTINUE 
CALL ERASE 
CALL POSAC(O) 

WRITE(6 , 100) 

100 FORMAT ( "ENTER CLASS NO."/) 

REA0(6,«) NUM 

READ IN SIGNAL IN CLASS SPECIFIED 

CALL EXECI 14,102B,ISG,128,NCLS,NUM) 

IDN0= ICD( NUM) 

FIND CLASS ENTRY IN MASTF AND READ INTO JBUF 

00 70 JM=l,Nv>'D 

IF ( JMFD( JM) +IDNO) 70,71,70 

70 CONTINUE 
GO TO 1 

71 CONTINUE 

CALL SL2RSI JMFO, 1 , JM, JR) 

L0CLI=JR 

CALL DSK2( 14, JBUF,384, LOCLI ) 

ZERO FILE MEANS AND VARIANCE 

DO 11 1=1,50 
JNPA( I ) =0 
CLMNI I ) =0. 

CLVA( I ) =0. 

11 CONTINUE 
NENT=0 

COMPUTE MEAN + VARIANCE OF PARAMETERS. 

DO 20 K=1 ,128 
IFdSGI KJ ) 21,21 ,22 

FIND RELATIVE SECTOR OF ISG 

22 DO 30 J=l,256 

IF( IPD( J J-ISG(K) ) 30,23,30 
30 CONTINUE 

WRITE(6,101) ISG(K) 

101 FORMAT ( "SIGNAL NO.", 16," NOT FOUND"/) 

GO TO 20 

23 NENT=NENT+l 

JGL0B( 3)=JGL0B( 3)+l 
NDX=2«J+2 

CALL EXECd4, 1 0 2B , S I G , 256 , NPAR , NDX ) 

CALL TRNSGd SGI K ) , IDNO) 

CALCULATE SUM + SUM OF SQUARES 

DO 40 1=1, 50 
A=SPAR( I ) 

GMI I ) = GMI I ) + A 
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GV( I } = GV( I )fA=i'A 
IF(JW(D) 40,40,50 

50 CONTINUE 

JNPA( I ) =JNPA( I } + I 
XNUM=FL0AT( JNPAI I i ) 

XNUM1=XNUM-1 . 

XNUM2=XNUM-2. 

CLMN( I ) = ( CLMNI I ) *XNUM1 +SP AR{ I ) ) /XNUM 
IFUNPAI I )-l) 51 ,51,52 
52 Continue 

CLVA( I ) =XNUM2/XNUM1«CLVA( I) + ( ( S P AR ( I ) -C LMN ( I ) )=!'=>=2 )/ 
1 XNUM 

51 CONTINUE 
40 CONTINUE 

20 CONTINUE 

21 CONTINUE 

CALL CURSRI IC, IX, lY) 

DO 80 1=1,50 
JM=I+14 

JPNUMI I ) = JMF0( JM) 

80 CONTINUE 

CALL DSK2I 15,JBUF,384,L0CLI) 

CALL ERASE 
CALL POSAC(O) 

WRITE! 6,102) IDNQ 

102 FORMAT ("CLASS NAME " ,16/ 

1" PAR MEAN VARIANCE frUSED"/) 

WRITE! 6,103) ! I,CLMN! I ) ,CLVA! I ) , JNPA! I) ,1=1,50) 

103 FORMAT(2!I6,2E12.4,I4,2X) ) 

CALL CURSR! IC, IX, lY) 

IF!IC'*121B) 1,999,1 

999 CONTINUE 

CALCULATE GLOBAL MEAN AND STD DEV 

CALL ERASE 
CALL POSAC(O) 

WRITE!6,800) 

800 FORMAT ! "GLOBAL VALUES"/, 

1"PAR GLOBAL MEAN STD DEV TY P DEV JBND") 

801 FORMAT !6X, I 4, 3! 2X,E12.4),1X,I4) 

C= JGLGB!3) 

DO 802 1=1,50 
A=GM! 1 ) 

B=GV! I ) 

A=A/C 

D=SQRT ! ! B-C4'A*A ) / ( C-1 . ) ) 

IF(D) 803,803,804 

803 D=l. 

804 GM! I )=A 
GV! I )=D 
GTV! I ) = .2<'D 
JBND! I)=2 

802 CONTINUE 

SET UP VALUES IN SECTOR 0 OF MASTER FILE 

JMFO!5) =-l 
JMF0!6) =5 

CALL EXEC! 15, 102B, JMFO, 128 , NLIB,0) 

PRINT OUT VALUES OF GLOBAL FILE 

WRITEI6 ,801 ) ! JMFO! H-14) ,GM! I ) ,GV! I ) ,GTV ! I ) , 

1 JBND! I ) ,1=1,50) 

CALL DSK2! 15,JGLOB,384,LG) 

STOP 

END 
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SLBROUTINE TRNSG ( I S , LDUM) 

DIMENSION WSC(15) 

DIMENSION SPAR( 50) , JW( 50) 

DIMENSION NPAR(3) ,NCLS(3) ,NLIB(3) 

COMMON JSAVE(50) 

COMMON JMFO (128) ,JMFD(768), MFE(384) 

COMMON IP0(256) ,SIG( 128) 

COMMON ICO ( 128) , I SG( 128) 

COMMON JBUF(4) ,CLMW(50) ,CLVA(50) , JNPA(50) , JPNUM(50), 
1 JCD0(80) 

COMMON JGL0B(4) , GM(50) ,GV(50) ,GTV( 50) , JBND( 50) , 

1 JGDMI30) 

EQUIVALENCE ( SP A R ( 1 ) , MFE ( 5 ) ) , ( J W ( 1 ) , MFE ( 1 0 5 ) ) 

EQUIVAL ENCE ( I DNO, JBU F ( 1 ) ) , ( NEN T , J BUF( 3 ) ) 

THIS SUBROUTINE CALCULATES PARAMETERS AND STORES THEM 
IN THE MASTER FILE BUFFER MFE 



DATA WSC( 1) ,WSC( 2) »WSC(3) ,WSC(4) 

1 /lOOO. , 1 1000. , 1000./ 

DATA WSC(5) ,WSC(6) , H SC ( 7 ) , W SC ( 8 ) ,WSC(9) 

1 /I. , 100. , 1.E05, 1.E09, 1.E13/ 

DATA WSC( 10) ,WSC ( 1 1 ) , WSC ( 1 2 ) , WSC ( 13) . WSC (14) 
1 /I. , 100. , l.c06, l.E09,l .E13/ 



THIS SUBROUTINE SCALES PARAMETERS AND STORES THEM 
IN THE MASTER FILE BUFFER MFE 



NTCP=JMF0(2) 

DO 5 I=1,NT0P 
IF (IS-JMFD(D) 5,6,5 

5 CONTINUE 

DO 7 1=1,50 
JW( I)=0 
7 CONTINUE 
WRITE (6, 100) 

100 FORMAT! 15, IX, "NOT IN MASTF") 

RETURN 

6 CONTINUE 
K = 0 

CALL SL2RS( JMF0,1 , I, JR) 

CALL DSK2( 14 ,MFE ,384, JR) 

MFE(1)=IS 
MFE(3)=-JBUF( 1) 

MFE(4) =0 
K = K+1 
JWT = 1 

RATE 

SPAR(1)= SIG(1)*WSC(1) 

RESOLUTION 
SPAR(2)= SIG(2) 

POLYNOMIAL MEAN SQUARE FIT COEFFICIENTS 

SPAR(3)= SIG(5)*WSC(5) 
SPAR(4)=SIG(6)-WSC(6) 

SPAR(5)= SIG( 7)*WSC(7) 

SPAR(6)= SIG( 10)>:'WSC( 10) 

SPAR( 7) =SI G( 1 1 ) - WSC ( 1 1 ) 

SPAR(8)= SIG ( 12 ) =«=WSC( 12 ) 

LAGUERRE POLYNOMIAL COEFFICIENTS 
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SPAR(9)= SIG(15) 

SPAR(10)= SIG(16) 

SPAR(11J= SIG(17) 

SPAR( 12 )=S IG( 18 ) 

SPAR(13)= SIGU9) 

MEAN BIAS 

SPAR(14}= SIG(20) 

VARIANCE 

SPAR(15)= SIG(35) 

INTRINSIC VARIANCE 

SPAR(16)= SIG(36) 

NONPARAMETRIC CORRELATION 

SPAR(17)= 1.-SIG(36)/SIG( 35) 

FOURIER POWER SPECTRUM USING LINEAR INTERPOLATION 

K= 18 
J=40 
L = 1 

RES=SIG( 2) 

SPAR(K)=ABS(SIG{ J ) ) 

DO 60 I=l»32 
K=K+1 

63 CONTINUE 

DIFF = FL0AT( I ) *2 . -FLOAT ( L) * RES 
IF(OIFF) 61,61,o2 

61 SPAR(K) =SIG( J+1 )+( SI G(J+1)-SIG(J))«D IFF/RES 
GO TO 60 

62 L=L+1 
J = J + 1 

GO TO 65 
60 CONTINUE 
256 CONTINUE 

DO 70 I=1t50 
JW( I ) =1 
70 CONTINUE 

CALL DSK2( 15,MFE,384, JR) 

RETURN 

END 
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